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FOREWORD 


This  report  describes  the  method  and  computer  program  for  calculat¬ 
ing  laminar  boundary-layer  properties  over  bodies  of  revolution  at 
large  incidence  in  subsonic  flow.  With  minor  changes,  it  was  also 
Mr.  Lee  A.  Kama's  thesis  for  the  Master  of  Science  Degree  in  Mechanical 
Engineering  at  North  Carolina  State  University  in  1983. 
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Procurement  Instrument  Identification  Number  (Contract  Number)  F3361 5-81 - 
K-3625  with  the  Air  Force  Wright  Aeronautical  Laboratories,  Wright- 
Patterson  Air  Force  Base,  Ohio  45433-6553.  The  subject  contract  was 
initiated  under  Air  Force  Flight  Dynamics  Laboratory  Project  2307, 

Task  2307N324,  on  17  August  1981  and  was  effectively  concluded  in 
May  1983.  Mr.  William  H.  Lane,  AFWAL/FIGC,  was  the  Air  Force  Project 
Engineer  for  the  study.  Comments  may  be  directed  to  him  at  (513)255- 
8486,  or  in  writing  at  the  above  address. 
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SECTION  1 


INTRODUCTION 

There  are  at  present  numerous  methods  available  for  calculating 
boundary  layers  over  two-dimensional  and  axi symmetric  bodies  at  zero 
angle  of  attack  (Refs.  1  and  2).  Methods  for  calculating  three- 
dimensional  boundary  layers  (Refs.  3,  4  and  5)  are  not  as  numerous  due 
to  the  fact  that  they  have  only  recently  come  under  investigation. 

These  techniques  are,  at  present,  limited  in  application.  Generally, 
fully  three-dimensional  methods  require  considerable  storage  and  com¬ 
putational  time  on  existing  digital  computers.  To  compound  the  problem, 
the  potential  solution  in  most  cases  cannot  be  described  through  a 
simple  analytical  expression. 

A  relatively  simple,  approximate  method  for  calculating  three- 
dimensional  boundary  layer  properties  is  the  axi symmetric  analog.  In 
this  method  the  boundary-layer  equations  are  written  in  a  streamline 
coordinate  system  and  the  cross  flow  velocity  is  assumed  to  be  zero. 

This  reduces  the  three-dimensional  boundary  layer  equations  to  a  form 
that  is  identical  to  those  for  axi symmetric  flow,  provided  that  (1)  the 
distance  along  an  inviscld  surface  streamline  is  interpreted  as  dis¬ 
tance  along  an  "equivalent  axlsymmetric  body,"  and  (2)  the  metric  co¬ 
efficient  that  describes  the  spreading  of  the  streamlines  is  interpreted 
as  the  radius  of  the  equivalent  body.  This  allows  any  existing  axisym- 
metric  boundary- layer  program  to  be  used  to  compute  the  approximate 
three-dimensional  boundary-layer  properties.  By  considering  multiple 
streamline  paths,  an  entire  surface  can  be  covered. 


The  major  difficulties  in  applying  the  ax i symmetric  analogue  are 
the  calculation  of  the  inviscid  surface  streamlines  and  the  corresnonding 
scale  factor.  References  (7)  and  (8)  provide  two  methods  with  which  to 
trace  inviscid  surface  streamlines  from  surface  pressures.  These  two 
approaches  are  basically  identical  in  that  each  requires  the  integra¬ 
tion  of  a  first-order  differential  equation  to  yield  the  streamline  angle 
(the  angle  between  a  streamline  and  body  meridian). 

In  Reference  (7)  Vollmers  proposes  that  a  shooting  technique  be 
employed  to  determine  an  initial  value  for  the  streamline  angle.  He 
contends  that  a  valid  initial  value  may  be  obtained  near  the  stagnation 
point  if  the  streamline  angle  approaches  the  correct  limit  durinq  the 
upstream  integration  from  a  given  point  on  the  body.  For  a  nonspherical 
nose,  the  streamline  geometry  is  such  that  the  streamline  angle  is  either 
0  or  180  degrees  in  the  limit  at  the  stagnation  point  (Ref.  9).  Vollmers, 
however,  fails  to  account  for  the  behavior  of  the  scale  factor  during 
the  upstream  integration.  There  is  a  possibility  that  the  scale  factor 
may  tend  to  zero  during  the  integration  despite  the  fact  that  the  stream¬ 
line  angle  may  approach  the  proper  limiting  value.  A  scale  factor  of 
zero  implies  that  the  streamlines  cross  at  a  particular  point  and  this 
is  a  physical  impossibility.  This  condition  can  occur  when  approximate 
surface  pressures  are  used.  If  the  inviscid  surface  velocity  components 
were  known,  the  streamlines  could  be  calculated  more  easily  and  there 
would  be  no  possibility  of  streamline  crossing. 

In  Reference  (8),  DeJarnette  attacks  this  problem  by  describing  the 
streamline  geometry  analytically  in  the  stagnation  region.  Outside  this 
region,  the  inviscid  surface  streamlines  are  calculated  from  the  surface 
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pressure  distribution.  This  is  the  approach  employed  in  the  present 
method  for  instances  in  which  an  analytical  potential  solution  is  not 
available. 

The  method  used  to  calculate  the  "equivalent  radius"  along  a  stream¬ 
line  follows  from  a  method  developed  by  DeOarnette  in  Reference  (9)  and 
is  included  in  this  study.  DeJarnette  has  previously  used  a  two  stream¬ 
line  approach  to  determine  the  scale  factor,  but  this  necessitates  the 
calculation  of  a  second  or  auxiliary  streamline.  In  the  present  method, 
the  scale  factor  is  calculated  along  a  single  streamline  from  the  solu¬ 
tion  of  two  first-order,  auxiliary  differential  equations  which  are 
functions  of  the  surface  pressure  distribution  and  body  geometry. 

The  present  method  employs  Hall's  (Ref.  4)  and  Blottner's 
(Ref.  2)  methods  to  obtain  a  solution  to  the  axi symmetric  boundary-layer 
equations.  Hall  applies  a  Crank-Nicholson  differencing  technique  to  the 
nondimensionalized  equations.  The  body  radius,  which  appears  only  in 
the  continuity  equation,  is  replaced  by  the  scale  factor  in  accordance 
with  the  axisymmetric  analogue.  Blottner  employs  the  same  differencing 
technique  to  the  axisymmetric  boundary- layer  equations  written  in 
transformed  variables.  In  this  case  the  body  radius  appears  only  in 
the  definition  of  the  transformed  variables  and  is  likewise  replaced 
with  the  scale  factor.  The  body  radius  does  not  appear  explicitly  in 
the  transformed  boundary- layer  equations.  The  development  of  each  of 
these  methods  is  included  in  this  study.  In  the  development  of  the 
boundary-layer  code,  various  velocity  profile  convergence  tests  and 
boundary- layer  edge  tests  are  investigated  also. 
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Any  axisymmetric  configuration  may  be  input  to  the  program  as  long 
as  an  analytical  expression  for  the  pressure  distribution  is  provided. 

In  the  event  that  only  experimental  pressure  data  are  available,  the  input 
geometry  is  restricted  to  spherically  capped  bodies  due  to  the  limita¬ 
tions  of  the  techniques  used  to  represent  the  surface  pressure  distribu¬ 
tion.  On  the  spherical  cap  the  pressure  distribution  is  represented  by 
a  Fourier  cosine  series  while  on  the  remainder  of  the  body,  a  doubly 
quadratic  spline  is  used  to  model  the  experimental  pressure  data.  The 
body  geometry  may  be  expressed  in  either  English  or  SI  units,  or  in  non- 
dimensional  form. 

Results  from  the  computer  program  are  presented  for  a  sphere, 
ellipsoid  of  revolution  with  thickness  ratio  1/4  and  a  sphere-ogive- 
cylinder  configuration  as  example  applications  of  the  computational 
method. 
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SECTION  2 

INVISCID  SURFACE  STREAMLINES 

The  axi symmetric  analogue  concept  employed  in  the  computer  Drogram 
effectively  reduces  the  three-dimensional  nature  of  the  boundary  layer 
to  that  of  an  axisynmetric  one  along  an  inviscid  surface  streamline. 

This  necessitates  the  calculation  of  the  equivalent  radius  of  the  newly 
defined  axisymmetric  body.  The  equivalent  radius  or  scale  factor  is  the 
metric  for  the  coordinate  $  normal  to  the  streamline  on  the  body  surface 
and  is  calculated  along  inviscid  surface  streamlines.  The  scale  factor 
is  an  indicator  of  the  physical  spacing  between  streamlines.  A  scale 
factor  that  is  increasing  indicates  that  the  streamlines  are  diverging 
and  thus  the  equivalent  radius  is  increasing. 

The  method  of  DeOarnette  (Ref.  8)  is  used  to  trace  the  inviscid 
surface  streamlines.  In  this  method,  the  body  geometry  is  exoressed 
in  terms  of  the  unit  vectors,  @x,  @r  and  8^  which  form  an  orthogonal 
cylindrical  coordinate  system.  Unit  vector  §x  is  parallel  to  the  body 
axis,  unit  vector  §r  is  in  the  radial  direction  and  normal  to  the  body 
axis.  The  third  unit  vector  in  this  system,  8^,  is  in  the  circumferential 
direction  (see  Figure  2.1). 

A  second  coordinate  system  which  is  oriented  to  the  body  surface 

is  used  to  describe  the  surface  streamlines.  This  system  consists  of 

the  unit  vectors,  §,,,  8  ,  and  §,.  Unit  vector  8  is  normal  to  the  body 

li  n  $  n 

surface  and  is  given  by 

8  =  -  sin  r  8  +  cos  T  8„  . 
n  x  r 
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(2.1) 
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Unit  vector,  e^,  is  tangent  to  the  body  surface  and  lies  in  a  meridianal 
plane.  This  vector  may  be  expressed  as 

gjj  =  cos  r  ix  +  sin  r  §r  .  (2.2) 

The  angle  r  is  the  body  angle  and  is  a  function  of  x  only  for  an  axisym- 
metric  body  (see  Figure  2.2). 

The  streamline  geometry  is  expressed  in  terms  of  unit  vectors, 
e^  and  en  which  also  form  an  orthogonal  coordinate  system.  Unit  vector 
§s  is  along  a  streamline  and  tangent  to  the  body  surface.  Unit  vector 
gg  is  normal  to  the  streamline  and  also  tangent  to  the  body  surface. 

Unit  vector  @n  is  used  in  common  with  the  previous  coordinate  system. 
Since  the  streamline  is  projected  on  the  body  surface,  the  component 
along  §n  is  zero  and  a  streamline  will  thus  lie  in  the  tangent  plane 
generated  by  unit  vectors  @„  and  §fl.  This  is  the  identical  plane  qener- 
ated  by  e^  and  g^.  The  streamline  angle  9  is  defined  to  be  the  angle 
between  unit  vectors  e$  and  g^  (see  Figure  2.3).  This  angle  is  the 
inclination  of  the  streamline  relative  to  a  body  meridian.  /The  stream¬ 
line  unit  vectors  may  then  be  written  in  terms  of  the  body  geometry  unit 
vectors  as 

@s  =  (cos  0  cos  I)  @x  +  (cos  0  sin  r)  @r  +  (sin  0)  e^  (2.3) 

and 


gg  =  -(sin  0  cos  r)  gx  -  (sin  0  sin  r)  gf  +  (cos  0)  g^  .  (2.4) 
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Figure  2.3.  Streamline  Coordinate  System 
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DeJarnette  (Ref.  9)  then  constructs  the  transformation  operators 
which  relate  the  streamline  partial  derivatives  to  the  cylindrical  co¬ 
ordinate  derivatives.  The  operators  may  be  written  as 


1  J_ 
h  98 


(e„ 


1A 


JL 

30 


(2.5) 


and 


JL 

3d> 


(2.6) 


where  D/DS  is  a  derivative  along  a  streamline  and  ^  ~  is  a  derivative 

n  9p 

normal  to  a  streamline  and  on  the  body  surface.  Substituting  the 
expressions  for  the  unit  vectors  yields 


and 


I  JL 

h  38 


sin  6  cos  r 


_3_ 

3x 


cos  e  jl 
r  9<j> 


(2.7) 


Application  of  Equation  (2.8)  yields 


Dx 

DS 


COS  0  cos  r 


(2.8) 


(2.9) 


and 


0*  _  sin  9 
DS  ‘  r 


(2.10) 
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These  differential  equations  may  be  numerically  integrated  to  give  the 
axial  and  circumferential  position  along  a  streamline  when  6  is  known 
(the  initial  values  for  all  differential  equations  will  be  discussed 
later  in  this  section).  If  the  potential  solution  were  known  in  analytic 
form,  the  angle  0  could  be  determined  from  it.  When  only  pressure  data 
are  available,  6  must  be  calculated  from  a  differential  equation. 

DeJarnette  derives  the  differential  equation  for  the  streamline 
angle  0  from  the  application  of  Euler's  equation  on  the  surface  of  the 
body.  In  vector  form  Euler's  equation  is 


DU  _  Vp. 
Dt  “  "  p 


(2.11) 


The  convective  term  may  be  recast  in  streamline  coordinates  as 


D§ 

ubK  +  u!-d#  •  <2-12> 

The  pressure  gradient  in  streamline  coordinates  may  also  be  written 
as 


VD  =  g  +  I  £  +  lE  @ 
vp  DS  es  h  30  e0  3n  en 


(2.13) 


Euler's  equation  may  then  be  written  as 


U  S  §  +  U 


Q* 

2  _ s  _  i/D£fl  +1l£a  +1£a\ 

DS  cs  '  u  DS  '  "  p  \DS  @s  +  h  36  80  +  3?  §n/ 

The  scalar  product  of  e^  with  this  equation  gives 


(2.14) 


D§e  , 

n2  — —  .  §  =  _  -1.  JlB. 
u  DS  e0  ph  30  * 


(2.15) 
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T’ ?  scalar  product  of  with  the  derivative  of  the  expression  for  e$ 
yields  (Ref.  9) 

d£  •  §  =  d0  +  sin  r  d<| ;  .  (2.16) 

S  p 

Substitution  of  this  equation  into  (2.15)  yields 


D0 

DS 


&  +  sin  r 


Dcp  _  _ 1  1^  9j) 

OS  "  '  pU‘  h  96  ■ 


(2.17) 


Application  of  the  previously  defined  transformation  operators  can  be 
used  to  write  this  equation  as 


06  _ 
DS  " 


1 

2 


■N 

u 

2 

f  9C  _ „ 

9C 

GO 

lUej 

« 

sin  0  cos  T  3x  +  r 

D 

8(f) 

4 

sin  r  sin  0 


(2.18) 


Equations  (2.9),  (2.10)  and  (2.18)  can  be  numerically  integrated  to  fully 
describe  the  geometry  of  the  streamlines  resulting  from  a  given  pressure 
distribution.  If  the  potential  solution  is  expressed  analytically,  0 
may  also  be  expressed  analytically  and  only  Equations  (2.9)  and  (2.10) 
need  be  integrated. 

The  scale  factor  must  be  evaluated  simultaneously  during  the  stream¬ 
line  integration  for  use  in  the  boundary- layer  calculations.  The  tech¬ 
nique  employed  here  has  not  previously  been  used  elsewhere  but  follows 
from  a  technique  developed  by  DeJarnette  (Ref.  9).  If  the  streamline 
coordinate  8  is  substituted  into  the  transformation  operators  (2.7)  and 
(2.8),  the  result  is  easily  shown  to  be 


and 


sin  0  cos  r  || 


cos  6  98 
r  9<f> 


(2.19) 
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COS  0  COS  r  || 


Dfi 

respectively.  Note  that  nf  =  0  since  the  coordinate  0  is  constant  along 


+  sin  e  33 
r  30 


(2.20) 


DS 

a  streamline.  Equation  (2.20)  may  be  solved  for  — 
into  Equation  (2.19)  to  yield 

sin  0  cos  r 


and  then  substituted 


I 

h 


-sin  0  33 

1  .  COS  0  33 

r  cos  0  cos  r  30 

XJ 

r  3* 

(2.21) 


which  reduces  to 


I  _  1  33 

h  “  r  cos  0  30 


(2.22) 


or 


h  = 


r  cos  0 


(2.23) 


Note  that  since  3  3  B(x,0), 


33 

30 


x 


This  expression  necessitates  the  additional  calculation  of  If- 

OP 

a  streamline  in  order  to  calculate  the  scale  factor  h. 


along 


The  differential  equation  for 


30 

33 


may  be  obtained  as  follows. 


Equations  (2.9)  and  (2.10)  may  be  combined  to  give 


D£  a  tan  0 
Ox  ”  r  cos  r  ' 


(2.24) 
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Taking  ^ 


of  both  sides  of  this  equation  yields 


J_  M  u  sec2  9  30 

3 <p  (DxJ  “  r  cos  T  S<p 


(2.25) 


since  both  r  and  r  are  functions  of  x  only  for  an  axisymmetric  body. 
With  8  =  8(x,<j>), 


S_ 

dip 


=  !£ 

3 

x 

38 

x 

(2.26) 


and 


D4  _  _3£ 
Dx  '  3x 


(2.27) 
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the  derivative  of  Equation  (2.24)  with  respect  to  <p  may  be  written  as 

(2.28) 


3 

34 


This  in  turn  may  be  rewritten  as 


fD4l 

,  38 

3 

_  sec2  e  30 

(DxJ 

3<)) 

x  36 

X 

3x; 

g  r  cos  r  3(j> 

_3_ 

fPi]  ,  _J 

D 

[34] 

3<t> 

[DxJ  3£ 

Dx 

36j 

38 


(2.29) 


or 


_D_ 

Dx 


fin  || 

■ 

38 

1  30 

r  COS2  8  COS  r  dip 


(2.30) 


By  application  of  the  chain  rule  and  Equation  (2.9),  this  differential 
equation  may  be  recast  as 
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D 

DS 


in 

i  ' 

1  38 

w  38 

'x 

r  cos  0  dcp 

(2.31) 


The  integration  of  this  differential  equation  along  a  streamline  requires 
30  I 

that  -rr  also  be  known. 

3<?>L 

*  30 

The  differential  equation  for  may  be  derived  as  follows:  sub- 

^  x 

stitution  of  the  transformation  operator  In  Equation  (2.8)  into  the 
streamline  Equation  (2.18)  yields 


rn«.  rn«.  r  39  +  sin  6  30  _  sin  8  sin  r  _  1 
cos  6  cos  r  3x  +  r  -  r  ? 


u 

2 

[15>1 

00 

h  38  ^ 

(2.32) 


30 


Solving  for  ^  gives 


30  _ _ 1 

3x  cos  0  cos  r 

Taking  the  partial  derivative  with  respect  to  <j>  of  Equation  (2.32) 
yields 


sin  8  30  sin  T  sin  8  1^ 

r  34>  r  2 


u 

2 

iM 

00 

"ue 

[h  36  J 

A 

(2.33) 


cos  8  cos  r  ^30x  +  S1"  0  |p-  -  sin  0  cosT  ||  ||  + 


30 

(34>. 


sin  T  cos  8  30  t  1 
r  H  2 


u 

2 

00 

u 

ej 

coseco5rS|i.il^^|i 


+  sin  0  cos  r 


cos  8  3Cp^ 
~  r  3<p 


32C. 


3<£3x 


a2c 

cos  8  P 

fu 

oo 

2  J_!!!e 

r  3<r 

Ufi  d<p 

3C 

sin  8  cos  r  -jf. 


(2.34) 
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Bernoulli's  equation  enables  the  circumferential  velocity  derivative  in 
this  expression  to  be  written  as 


ue  3$ 


1 5>  ty 

‘  2  34)  u  I 


(2.35) 


Substitution  of  this  equation  and  the  two  transformation  operators  into 
Equation  (2.34)  gives 


D  (901  .  n  rrir  r  90  90  COS  8 

DS  3?  =  sin  6  cos  r  ta  H  -  ~r~ 


1  fOa  DC  9*C 

+  2  ~DS  30  +  Sin  8  C0S  f  9*9^ 


,  fuj"  9C  f.  3C 

I  _ g.  II _ g> 


2  ue  34)  (h  30 


sin  r  cos  8  90 
r  3<{> 


COS  9  9  Cp 
r 


(2.36) 


Substitution  of  Equation  (2.33)-  into  this  equation  yields 

D  f 90  _  sin  8  [39  sin  8  [38  sin  r  sin  9 

DS  3cj>  "  "  cos  0  3(j>  r  3q>  r 

XJ  XJ  _  xd 


+  JZ  _1_  cos  8  [30 

ue  2h  98  r  3<j> 


38 1  l1 2  _  sin  r  cos  0  [38 


,  1  N2  0C-  [39  C„-  r  82C»  cos  e  3!cp 

*2  r  k  i*.  ti,n0msrw;--rF 
e)  x) 


yj  _  v  _ 

94>9x  r  34> 


1  fuJ4  3C  [.  3C 

1  —  P  I  P 

2  ue  94>  (h  98. 


(2.37) 
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Note  that 


DC  •  a 
p  sin  p 

DS  '  cos  0 


15,' 

h  36 


3C  .  „  3C  -2a  gC 

cos  9  cos  r  J  *  5inl  ♦  HalS-SaiX  la 
3x  r  34>  cos  0  9x 


sin  8  3^p 
r  3<J> 


This  may  be  simplified  and  written  as 


0Cp  _  sin  e 
DS  "  cos  0 


i5>l 

h  36 


_  cos  r  9Cp 
cos  e  3x 


Substitution  of  this  expression  into  Equation  (2.37)  leaves 


D 

DS 


[30] 

sin2  6 

[30^ 

2 

sin2  6  sin  r 

fie] 

cos  0 

[30]  2 

W 

r  cos  0 

r  cos  0 

r 

M 

.  12  ap 

sin  r  cos  8  [30]  +  \  _»  [36]  cos  r  6  p  +  1_ 

r  [34>J  2  u  (30 J  COS  6  3x  2 


sin  0  cos  f 


3ZC  q  32C 

_ 2.  _  cos  9  D 

34>3x  ’  r  3<fr 


2 


y*^a 

ueJ  ^ 


'u  ]2 
00 


i5>] 

h  36  „ 


(2.38) 


With  the  substitution  of  the  transformation  operator  in  (2.7),  Equation 
(2.38)  may  be  rewritten  as 


D_ 

DS 


[39 

3  ip 

[30 

30 

+  sin  r| 

3  <p 

x) 

<t<P 

x _ / 

r  cos  0 


*  i  N 


N  L 


sin  0  cos  r 


3<f>3x 


a  32C 

COS  0  p 

+  I 

[30 

uJ 

00 

1 2  3C 

1  cos  r  %  ,  l 

Uco 

r  34» 

+  2 

dt 

X 

ueJ 

j  cos  0  3x  2 

Lue 

4  5, 

3<f> 


sic  0  cos  r  ^ 

r  3<p 


3x 


(2.39) 
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This  differential  equation  involves  first  and  second  derivatives  of  the 
pressure  coefficient  which  are  supplied  by  either  the  analytical  solu¬ 
tion  or  the  spline  fit  if  only  experimental  pressures  are  supplied. 

a2c 

Note,  however,  that  —  - j  does  not  appear  and  thus  the  spline  fit  in  the 

oX 

x-direction  is  simplified.  Equations  (2.31)  and  (2.39)  may  be  integrated 
to  give  the  scale  factor  (Equation  (2.23))  along  an  inviscid  surface 
streamline.  The  differential  equations  are  singular  at  the  stagnation 
point;  therefore,  the  geometric  position  where  the  integration  of  the 
streamline  equations  begins  must  be  some  distance  away  from  the  stag¬ 
nation  point  (Ref.  9).  The  technique  used  here  depends  on  whether  the 
inviscid  properties  are  obtained  from  experimental  pressure  data  or 
an  analytical  potential  solution. 

When  experimental  pressure  data  are  used,  there  are  generally 
insufficient  data  near  the  stagnation  point  to  adequately  determine 
the  pressure  distribution  needed  for  the  Integration  of  the  streamline 
equations  near  the  nose.  A  potential  panel  method,  USSAERO,  was  used 
to  obtain  additional  pressure  data  in  the  nose  region.  It  was  found 
that  this  pressure  data  and  the  experimental  pressure  data,  at  the 
most  forward  position,  were  reasonably  close  to  a  spherical  pressure 
distribution  about  the  stagnation  point.  A  spherical  pressure  distri¬ 
bution  produces  streamlines  along  spherical  meridians  from  the  stagna¬ 
tion  point  to  the  sphere-afterbody  Interface.  (See  Figure  4.1  on  page 
35.)  On  the  spherical  cap,  the  streamline  geometry  and  metric  are  given 
in  Ref.  74.  Integration  of  the  streamline  differential  equations  begins 
at  the  sphere-afterbody  Interface.  For  a  given  circumferential  angle  <f>  , 
page  31  of  Ref.  9  gives  the  initial  streamline  slope  as 
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eo  =  cos-1  |(cos  <*eff  cos  Tj  -  sin  ag^  sin  Tj  cos  <j>0)/sin  ^0|  (2.40) 

where  the  spherical  angle  is  determined  from 

ipo  =  cos-1  jcos  sin  Tj  +  sin  cos  Tj  cos  <f>Qj  .  .  (2.41) 

The  initial  value  of  (30/3<|>)o  for  Equation  (2.39)  is  obtained  by 
differentiating  Equation  (2.40).  The  result  is 


90_ 


o 


sin  a  ~~  sin  r.  sin  <j>  -  cos  0  cos  tpn 
_ eff _ I _ _o _ o _ _o 

sin  tpQ  sin  0Q 


(2.42) 


From  Ref.  14,  the  streamline  metric  on  the  spherical  cap  is  given  by 
h  =  Rper  sin  ip;  and  since  r  =  Rpgr  cos  T,  Equation  (2.23)  gives 


dip 

30 


sin 

_ To 

COS  \  COS  0„ 
0  0 


(2.43) 


as  the  initial  value  for  Equation  (2.31)  at  the  sphere-afterbody  inter¬ 
face. 

When  an  analytical  potential  solution  Is  known,  the  inviscid  ve¬ 
locity  components  can  be  used  to  obtain  an  analytical  expression  for 
the  streamline  angle,  0,  and  Its  circumferential  derivative,  30/3x. 

Then  Equation  (2.39)  is  not  needed  and  Equations  (2.9),  (2.10),  and 
(2.31)  can  be  integrated  numerically  to  determine  the  streamline  loca¬ 
tion  and  metric.  Initial  conditions  for  the  streamline  location  are 
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determined  from  an  axial  position  and  circumferential  angle  near  the 
stagnation  point.  As  discussed  in  Ref.  9,  the  initial  value  of 
( 34>/33) Q »  and  hence  hQ,  for  Equation  (2.31)  is  arbitrary.  Since  Equa¬ 
tion  (2.31)  is  used  to  integrate  2.n(3<J>/aB)0  will  have  no  effect  on  the 
numerical  integration  of  this  differential  equation.  The  actual 
value  calculated  for  3<J>/33 ,  however,  will  be  relative  to  the  initial 
value  (9<j)/3$)0. 


SECTION  3 


BOUNDARY  LAYER  METHODS 


The  present  computer  program  has  the  option  of  employing  either 
Blottner's  (Ref.  2)  or  Hall's  (Ref.  4)  boundary  layer  method.  For  each 
method  a  solution  is  obtained  through  the  use  of  a  finite-difference 
technique.  After  application  of  the  respective  transformation,  the 
governing  equations  are  then  cast  in  second-order  accurate  finite- 
difference  form.  Since  the  governing  equations  are  parabolic,  the 
boundary  layer  may  be  calculated  by  "marching"  downstream  in  a  step- 
by-step  fashion  along  an  inviscid  surface  streamline. 

Blottner's  method  involves  solving  the  governing  equations  written 
in  F-V  similarity  form.  These  equations  are  obtained  by  the  application 
of  the  Levy- Lees'  transformation  defined  for  incompressible  flow  as 


C(s)  3  Kpy  u 


r  S  U. 


„  UM 
0  00 


r2  ds 


and 


(3.1) 


n(s,n) 


ue  rpn 


rr  . 


(3.2) 


In  the  axisymmetric  analogue,  the  body  radius  r  is  replaced  by  the  scale 
factor  h  and  s  is  distance  along  the  inviscid  surface  streamline.  This 
transformation  creates  a  (E.n)  computational  grid  from  the  (s,n)  physi¬ 
cal  grid.  The  computational  grid  has  been  effectively  stretched  in 
both  the  normal  and  tangential  directions.  The  resulting  equations  are 
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(3.3) 


2Z 


0 


and 


2EFfftvf ♦  6(Fi  - 11 


(3.4) 


where 


V 


2Z 


F  |J  +  pvh//Tp</Kpuueh2 


(3.5) 


F  =  u/ue 


(3.6) 


and 


6  ue  d£ 


(3.7) 


(see  Appendix  A.l)  .  These  equations  are  then  cast  in  finite  difference 
form  using  the  Crank-Nicholson  scheme  to  yield  a  system  which  is  second- 
order  accurate  in  both  spatial  directions.  The  resulting  system  can 
be  conveniently  written  as 


A2  Fi+l,j-l  +  B2  Fi+l,j  +  C2  Vi+4,.j-l  +  E2  Vi+H,.j  =  °2 


where 


(3.8) 


A2  =  AnU  + 

B2  -  An(^  +  £i+>./AC) 
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and 


for  the  continuity  equation  and 


A1  F1+l,j-l  +  B1  Fi+l,j  +  C1  F1+l,j+l  +  E1  Vi+Js,j  =  D1  (3’9) 


where 


A,  ■  -  %  0  ♦  *n 

Bj  =  1  +  An2  ^+1  j  + 

ci  -  -  \  0  -  \  An  V1+Jj>j) 

E1  s  \  An  (Fi,j+l  "  Fi  »j-l  +  Fi+lJ+l 


Fi+l,j-l} 


and 


°1  ■  *  '  2F1,j  4  FlJ.J>  +  *"*  eUH  [(1  +  FW,j> 

4  (1  '  Fl,jO  *  ^  7f+%. J  Ijf+l.j-H  '  F1*l.j-l] 


+  An2  Cf+Js  (Fj+1j  +  F5j)/AE 


for  the  momentum  equation  (see  Appendix  A. 2).  To  provide  that  only 
a  linear  system  of  equations  needs  be  solved  to  obtain  F,  the  nonlinear 
terms  in  the  finite  difference  expressions  have  been  linearized  using 
the  Newton- Raphson  technique.  This  will  necessitate  repeated  iteration 
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in  order  to  achieve  a  converged  solution  to  the  actual  nonlinear  equa¬ 
tions.  The  bar  indicates  quantities  from  a  previous  iteration. 

By  virtue  of  the  transformation  in  Equations  (3.1)  and  (3.2), 
Blottner's  method  may  be  applied  at  the  stagnation  point  to  yield  a 
limiting  velocity  profile.  In  addition,  a  similarity  solution  is  oos- 
sible  since  £  *  0  at  that  point  (see  Appendix  A. 3).  The  boundary 
conditions  consist  of  the  edge  and  wall  conditions.  At  the  boundary- 
layer  edge  the  condition 


F(n  =  ne)  =  1 

is  applied  at  each  step  of  the  Integration.  The  value  of  the  normal  co¬ 
ordinate  at  the  boundary  layer  edge,  ne#  must  be  provided  initially  and 
must  be  large  enough  to  account  for  the  entire  boundary  layer  thickness. 
The  no-slip  condition  demands  that 

F(n  »  0)  »  0 

at  each  step.  The  pressure  gradient  parameter,  B,  is  related  to  the 
velocity  gradient  and  for  spherical  flow,  6,  becomes  1/2. 

The  system  of  equations  may  then  be  solved  to  yield  the  limit  of 
F  at  the  stagnation  point  through  use  of  the  modified  Davis  algorithm 
(Ref.  2).  This  algorithm  solves  a  coupled  system  of  equations.  The 
profile  slope  at  the  wall  may  be  expressed  by  a  second-order  accurate 
expression  which  in  turn  is  used  to  evaluate  (Cf  . 

At  points  away  from  the  stagnation  region  it  is  necessary  to  solve 
the  complete,  nonsimilar  system  of  equations  (Equations  (3.8)  and 
(3.9))  which  Involve  the  transformed  step  size,  A£,  along  a  streamline. 
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With  a  prescribed  step  size  along  a  streamline,  the  expression  for  £  may 
be  numerically  integrated  for  subsequent  use  in  the  boundary  layer  equa¬ 
tions.  The  calculation  of  the  pressure  gradient  parameter,  which  is 
defined  as 


.  2£^e 

ue  dC 


is  evaluated  at  the  midpoint  of  the  computational  interval  along  the 
body  surface.  On  the  nose  region  of  a  spherically  capped  body,  this 
term  reduces  to 


i+h 


2 

1  due 
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ha  dip 

per 

ej 

i  +h 


and  on  the  afterbody  it  becomes 


u  )2  f3u  nY 
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ej 


+ 


I!!*  m]| 
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The  derivative  and  both  and  are  supplied  by  subroutines 
SPHCAP  and  INVISO,  respectively.  The  two  total  derivatives  along  a 
streamline,  and  ,  are  used  In  conjunction  with  the  streamline 
integration  (see  INVISCID  SURFACE  STREAMLINES).  The  nonsimilar  equa¬ 
tions  may  then  be  solved  for  F  using  the  same  computational  technique 
as  was  used  at  the  stagnation  point. 

Hall's  method  involves  solving  the  governing  eauations  written  in 
terms  of  dimensionless  primitive  variables.  Hall  employs  the  customary 
transformation 
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s*  =  s/L 


"*  *  f 

r*  =  r/L 

u*  =  u/u 
'  00 

V*  =  v/U^ 

where 


which  yields 


3(u*r*) 

3s* 


+  r* 


for  the  continuity  equation  and 


(3.10) 


u* 


1  d(ue*) 

2  ds* 


2 


+ 


32U* 

an*^ 


(3.11) 


for  the  s-momentum  equation  (see  Appendix  A. 4).  The  desired  unknown 
is  u*  which  is  the  dimensional  velocity  normalized  by  the  frees tream 
velocity.  These  equations  are  then  cast  in  second-order  accurate  finite- 
difference  form  (see  Appendix  A. 5).  The  resulting  system  may  be  ex¬ 
pressed  as 


A2  Vl,j  +  B2  VlJ-1  +  C2  vi+>s,j  +  E2  Vi+*5,.i-l  =  D2  (3-12) 


where 
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and 


c2  =  1 


E2  =  -! 
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for  the  continuity  equation  and 


where 


A1  u1+l,j-l+  B1  ui+l,j  +  C1  u1+l,j+l  +  E1  V1+Js,j  s  D1 
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for  the  momentum  equation  (note  that  the  stars  have  been  omitted  for 
clarity).  These  equations  have  also  been  linearized  using  the  Newton- 
Raphson  method.  The  system  of  equations  which  results  at  each  station 
along  the  body  is  block  tridiagonal  in  form  and  may  be  easily  and  ef¬ 
ficiently  solved  in  the  same  manner  as  was  used  for  Blottner's  method. 

The  boundary  conditions  required  for  a  solution  to  Hall's  equa¬ 
tions  are  that 


u(s,nfi)  =  ue 


and 


u(s,o)  =  0 


for  the  no-slip  condition. 

To  begin  the  Integration  of  the  boundary  layer,  an  initial  profile 
must  be  known.  The  stagnation  point  is  an  ideal  place  to  start  the 
Integration.  Hall  (Ref.  4)  and  Geissler  (Ref.  10)  both  utilize  the 
well-known  three-dimensional  stagnation  point  boundary -layer  solution 
of  Howarth  (Ref.  11).  This  is  an  unwarranted  complication  since  at 
the  stagnation  point 

u(o,n)  =  0 
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at  all  points  across  the  boundary  layer.  This  may  be  used  as  the  initial 
velocity  profile  in  Hall's  method,  and  then  the  first  station  away  from 
the  stagnation  point,  along  an  inviscid  surface  streamline,  may  be  cal¬ 
culated  by  the  method  given  above. 


Convergence  Criteria 


Since  the  finite-difference  equations  have  been  linearized,  repeated 
iteration  is  necessary  in  order  to  obtain  a  solution  to  the  nonlinear 
equations.  The  iterative  process  could  be  made  to  continue  until  the 
solution  becomes  exact  (within  the  accuracy  of  the  computer)  but  this  is 
no  doubt  unwarranted.  In  practical  applications,  the  iterative  process 
is  usually  allowed  to  continue  until  the  solution  is  changing  by  less 
than  a  prescribed  amount  between  successive  iterations.  This  is  one 
definition  of  a  converged  solution. 

Because  the  skin  friction  Is  the  one  of  the  more  important  parameters 
of  interest,  it  appears  logical  that  convergence  should  be  based  on  it. 

In  practical  applications,  the  Iterative  process  should  stop  when  the 
skin  friction  changes  by  less  than  a  prescribed  amount  between  successive 
iterations.  This  is  the  definition  most  commonly  applied  in  two-dimensional 
boundary  layer  cases. 

The  computational  method  developed  has  the  option  of  employing 

either  of  these  definitions.  The  Input  parameter,  NC,  corresponds  to 

the  method  which  Is  used  to  define  a  converged  solution.  The  option 

corresponding  to  NC  =  0  specifies  that  convergence  Is  based  on  Cf  Re  J 

LI  e 

changing  by  less  than  0.5  percent  between  successive  Iterations. 

Covergence  is  based  on  the  velocity  at  each  grid  point  changing  by 
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less  than  0.1  percent  between  iterations  when  NC  =  1.  Table  3.1 
illustrates  the  effect  of  this  option  on  the  separation  point  for 
a  sphere  in  incompressible  flow. 

Table  3.1.  Effect  of  Convergence  Criteria  on  Boundary-Layer  Separation 
for  Hall's  and  Blottner's  Methods  on  Sphere  in  Incompres¬ 
sible  Flow 


Method 

NC 

Separation 

Angle  (Deg.) 

Steps  to 
Separation 

Hall 

0 

107.43 

184 

As  =  0.01 

An  =  0.0471 

1 

104.83 

184 

Blottner 

0 

105.75 

278 

AC  =  0.005 

An  =  0.11539 

1 

104.94 

276 

40  Points  Across  Boundary  Layer  Initially 


3.2.  Boundary  Layer  Edge  Criteria 

Since  the  velocity  in  the  boundary  layer  only  approaches  the  value 
of  the  inviscid  stream  asymptotically,  an  effective  edge  must  be  imposed. 
The  velocity  at  the  grid  point  which  is  arbitrarily  said  to  lie  at  the 
edge  is  assigned  the  velocity  of  the  inviscid  stream.  The  relationship 
of  the  velocities  at  the  grid  points  In  the  region  near  the  imposed 
edge  may  then  be  used  to  assess  whether  the  actual  boundary- layer  thick¬ 
ness  has  been  adequately  accounted  for.  According  to  the  classical 
definition,  the  boundary  layer  thickness  Is  adequately  reoresented  if 
the  velocity  at  the  point  adjacent  to  the  Imposed  edge  is  a  certain 
percentage  of  the  velocity  at  the  Imposed  edge.  This  percentage  is 


usually  In  the  range  of  99.5  to  99.995  percent.  Wang  (Ref.  5)  emoloys 

this  definition  in  his  fully  three-dimensional  technique.  A  second  test 

for  the  boundary- layer  edge  could  be  constructed  which  utilizes  the  fric- 

cf  yf^e,  ,  at  the  edge  of  the  boundary  layer  as  the 
k  e 

governing  criterion.  If  this  parameter  is  below  a  prescribed  limit,  the 
imposed  edge  may  be  considered  to  adequately  account  for  the  boundary- 
layer  thickness. 

Both  of  the  tests  described  above  are  Included  as  options  in  the 
present  computational  method.  The  parameter  NT  specifies  an  option  to 
be  used  for  the  edge  test.  The  option  NT  =  0  specifies  that  the  edge 
test  be  based  on  the  classical  definition  in  which  the  tolerance  is 


tion  parameter. 


99.95  percent.  The  second  test,  which  corresponds  to  NT  =  1,  requires 


hJv), 


be  less  than  0.005  for  the  boundary- layer  thickness  to 


be  adequate.  Table  3.2  illustrates  the  effect  of  both  the  edge  test 


and  convergence  test  options  on  the  separation  point  for  a  sphere  in 


incompressible  flow.  Note  that  the  more  stringent  option,  NT  =  1,  re¬ 


sults  in  the  addition  of  points  at  the  boundary-layer  edge  and  a  more 
accurate  separation  point.  The  results  generated  in  conjunction  with 
the  option  corresponding  to  NT  =  0  could  most  likely  be  improved  if  the 


respective  tolerance  were  to  be  decreased.  It  is  evident  from  the  re¬ 
sults  obtained  with  Hall's  method  that  the  boundary  layer  is  undoubtedly 
thickening.  Since  Hall  makes  use  of  primitive  variables,  a  growing 
boundary  layer  will  require  that  the  outer  edge  be  adjusted  occasionally. 
The  transformed  normal  coordinate  used  in  Blottner's  method  has  pro¬ 


visions  to  account  for  the  growth  of  the  boundary  layer.  Because  of 
this,  it  is  seldom  necessary  to  manually  shift  the  outer  edge. 
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Table  3.2.  Effect  of  Edge  Test  and  Convergence  Test  on  Separation 
for  Hall's  and  Blottner's  Methods  on  Sphere  in  Incom¬ 
pressible  Flow 


Method 

NC 

NT 

Separation 
Angle  (Deg.) 

Steps  to 
Separation 

Points 

Added 

Hall 

0 

0 

108.59 

194 

23 

As  =  0.01 

An  =  0.471055 

0 

1 

107.43 

184 

31 

1 

0 

104.89 

184 

18 

1 

1 

104.83 

184 

30 

Blottner 

0 

0 

105.02 

278 

0 

AC  =  0.005 

An  =  0.11539 

0 

1 

105.75 

278 

5 

1 

0 

104.94 

276 

0 

1 

1 

104.94 

276 
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SECTION  4 


SURFACE  PRESSURE  DISTRIBUTION 

If  an  analytical  potential  solution  Is  not  available  for  a  particular 
configuration,  experimental  pressure  data  must  be  applied.  The  accurate 
surface  fitting  of  the  pressure  data  Is  critical  not  only  to  the  calculation 
of  the  boundary- layer  properties  but  also  to  the  calculation  of  the  Invlscld 
surface  streamlines.  Near  the  nose  region  of  a  body  where  pressure  gradients 
are  relatively  large,  experimental  pressure  data  are  generally  not 
available.  The  region  of  the  body  downstream  of  the  nose  generally 
experiences  more  moderate  pressure  gradients  and  sufficient  experimental  data 
are  provided  to  model  a  surface  pressure  distribution.  After  Investigating 
several  methods  for  surface  fitting  experimental  pressure  data.  It  was  found 
that  a  doubly  quadratic  spline  would  adequately  model  the  pressure 
distribution  downstream  of  the  nose. 

As  mentioned  earlier,  a  potential  panel  method,  USSAERO,  was  used  to 
calculate  additional  pressure  data  In  the  nose  region.  Attempts  to  use  a 
doubly  quadratic  spline  to  blend  the  pressures  calculated  by  the  USSAERO  code 
In  the  nose  region  with  the  experimental  data  downstream  were  unsucessful. 

To  model  the  pressure  distribution  In  the  nose  region,  an  alternate  approach 
was  employed.  This  method  has  been  tested  for  a  sphere-ogive-cylinder  only. 
First,  the  pressure  data  calculated  by  the  USSAERO  code  were  blended  with  the 
experimental  data  at  the  most  forward  station  and  then  they  are  plotted  as  a 
function  of  the  spherical  angle  i  about  the  stagnation  point  (see  Figure  5.2 
on  page  42).  It  was  found  that  the  Fourier  cosine  series 
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(4.1) 


9 

Cn  M  =  An  +  l  An  C0S  NO 

p  0  n=1  n 

represented  the  pressure  distribution  in  the  nose  region  quite  satis¬ 
factorily.  In  this  series,  the  spherical  angle  <p  is  given  by 

ip  =  cos-1  jcos  sin  r  +  sin  cos  r  cos  <f>j  (4.2) 

where  a  ^  is  the  angle  between  the  body  axis  of  symmetry  and  the  line 
which  passes  through  the  stagnation  point  (see  Figure  4.1).  Note  that 
ip  =  0  corresponds  to  the  stagnation  point. 

The  coefficients  in  the  Fourier  series,  An  (n  =  0,9),  are  obtained 
from  the  solution  of  ten  simultaneous  equations  generated  from  the 
application  of  Equation  (4.1)  to  ten  distinct  points  on  a  curve  faired 
through  the  pressure  data  calculated  from  USSAERO.  One  point  must 
be  the  stagnation  point  itself  which  is  determined  by  interpolating 
data  from  USSAERO. 

Away  from  the  nose  region,  pressure  gradients  usually  become  smaller 
and  a  doubly  quadratic  spline  may  be  used  to  fit  the  experimental  pres¬ 
sure  data  downstream  of  the  interface.  In  order  to  describe  the  doubly 
quadratic  spline,  consider  the  singly  quadratic  spline  first.  An 
interval 

X]  <  x  <  xN 

is  divided  into  N  subintervals.  Each  interior  subinterval  rage  is 
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Figure  4.1.  Geometric  Parameters 


X„  ,  +  X„  X  +  X.  . 

— ~  i  x  i  2  "  1  ("  =  2,  N  -  1) 


and  the  subintervals  on  the  left  and  right  boundaries  range  from 


x,  <  x  < 


X1  +  x2 


1  -  -  2 


and 


XN-1  *  XN 
2 


<  x  <  x 


N 


9 


respectively.  The  dependent  variable  at  each  of  the  points  xp  is  de¬ 
signated  by  yn  (see  Figure  4.2). 

There  exist  N-l  midpoints  in  the  total  interval.  The  midpoints, 

denoted  by  Xn>  may  be  computed  by  the  above  relations.  Corresponding 

to  each  of  the  midpoints  Xn  there  is  a  yet  undetermined  dependent 

variable  Y.  Each  of  the  Y  's  is  determined  such  that  there  is  con- 
n  n 

tinulty  of  the  function  and  its  first  derivative  between  adjacent  sub¬ 
intervals.  The  first  derivative  must  be  specified  on  the  left  and 
right  boundaries  of  the  Interval.  This  will  yield  a  system  of  N-l 
linear  equations  for  Yn  (n=l,N-l).  This  system  may  be  expressed  as 


fc. 

(bl  +  cl*  Y1  +  dlV2  '  — +  ylbl  *  y2  ^C1  *  dl' 

Vn-1  *  (bn  *  cn)Vn  *  d„Vl  '  yn  <*n  *  V  +  yn*l  <dn  +  cn)(n=2'N-2> 

*L 

+  yM_i  (*m_i  +  bM_,)  +  yM(cw  ,) 


'N 


aN-lYN-2  +  ^bN-l  +  CN-1^  YN-1  “  ‘  2  '  'N-l  '°N-1  ’  "N-l'  ’  JN"*N-1 
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(4.3) 


Figure  4.2.  Quadratic  Spline  Parameters 


where 


bj  =  2/Ax2 

an  "  AW[xn(AVl  *  4x„>] 
bn  '  [2  +  Axn/4xntl]/[Axntl  +  V 
cn  =  [2  +  4xn+2/Axn+l^ixn+2  *  4xrt+l^ 

dn  -  WtAWAlW  +  4Vl» 

%-l  '  2/Axn 
and 


This  system  forms  a  tri diagonal  matrix  and  the  unknowns  may  be  obtained 
through  use  of  the  Davis  algorithm  (Ref.  2). 

A  second-order  polynomial  about  the  point  xp  which  may  be  written 
as 


(x  -  x)2 

y(x)  ■  yn  +  y„(x  -  xn)  +  yJJ - (4.4) 

is  applied  at  the  single  data  point  which  lies  within  each  subinterval. 

This  equation  contains  only  two  unknowns  since  yn  is  known  at  the  data 

point  x  .  On  a  given  interior  subinterval,  y'  and  y"  are  related  to 
n  n  n 

the  dependent  variables  Yn_j  and  Yn  (which  have  already  been  determined) 
and  can  be  expressed  as 
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(4.5) 

(4.6) 

Thus  to  determine  the  value  of  the  dependent  variable,  y,and  its  de¬ 
rivatives  at  any  position  on  the  total  interval,  all  that  need  be  done 
is  to  determine  in  which  subinterval  the  independent  variable  lies. 

The  corresponding  coefficients  in  the  quadratic  expression  (Equation 
(4.4))  may  be  generated  from  Equations  (4.5)  and  (4.6). 

The  extension  of  the  one-dimensional  quadratic  spline  to  two 
dimensions  is  a  relatively  simple  process  and  is  performed  as  follows. 
One-dimensional  quadratic  splines  y(x,<J>|c)  are  formed  for  specified  values 
of  4>k.  For  a  given  value  of  x,  y(x ,4>k)  and  (x ,d>k )  are  calculated  for 
each  <j>k-  These  values  are  then  fitted  by  a  quadratic  spline  in  the  <J>  direc 
tion.  These  splines  can  then  be  used  to  calculate  y(x,4>),  |^-  (x ,<f>) 
and  (x,<|>)  for  a  given  value  <J). 

The  quadratic  spline  yields  a  function  which  is  continuous  and  has 
a  continuous  first  derivative.  The  second  derivative  Is  continuous  (and 
constant)  over  each  subinterval,  but  Is  not  constrained  to  be  continuous 
at  the  junction  of  the  subintervals.  It  Is  possible  for  inflection 
points  to  occur  only  at  these  junctions.  Should  an  inflection  point  be 
desired  at  a  specific  location.  It  may  be  Included  simply  by  the  addi¬ 
tion  of  two  data  points  such  that  the  midpoint  of  this  Interval  becomes 
an  endpoint  of  a  subinterval. 


"  4Vi  *  4,» 
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SECTION  5 


STREAMLINES  ON  SPHERICALLY  CAPPED  GEOMETRIES 

Experimental  pressure  data  were  obtained  on  a  sphere-ogive-cyl inder 
at  a  =  45°  from  the  first  row  of  pressure  taps  to  the  base  region.  No 
experimental  pressure  data  were  obtained  on  the  spherical  cap.  To 
assist  in  modeling  the  pressure  distribution  on  the  spherical  cap,  the 
USSAERO  potential  code  was  used  to  calculate  pressure  data  and  these 
data  were  interpolated  to  locate  the  stagnation  point.  Due  to  the 
large  pressure  variation  over  the  nose  region,  the  doubly  quadratic 
spline  function  used  to  model  the  pressure  distribution  downstream  of 
the  sphere-ogive  interface  was  found  to  be  unsatisfactory  for  the  nose 
region.  An  alternate  approach  described  in  SECTION  4  was  to  graph 
the  calculated  pressure  data  from  USSAERO  on  the  spherical  cap  as  a 
function  of  the  angle  \p,  given  by  Equation  (4.2),  which  is  the  spherical 
angle  measured  about  an  axis  passing  through  the  stagnation  point  and 
the  center  of  the  sphere  (see  Figure  4.1).  The  results  are  given  on 
Figure  5.2  and  they  show  that  the  pressure  distribution  is  reasonably 
close  to  a  spherically  symmetric  one.  With  the  assumption  of  a 
spherically  symmetric  pressure  distribution,  the  streamlines  on  the 
spherical  cap  will  simply  follow  spherical  meridians  about  the  axis 
through  the  stagnation  point  and  the  center  of  the  sphere.  The  angle 
between  this  axis  and  the  body  axis  Is  in  Figure  4.1  which  is 
quite  different  from  the  actual  angle  of  attack,  a. 

For  spherical  flow,  the  pressure  distribution  along  one  meridiarr 
is  indistinguishable  from  another  and  the  boundary  layer  is  truly 
axisymmetric.  Note,  however,  that  the  magnitude  of  the  pressure  over 
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Figure  5.2.  Spherical  Symmetry  of  Nose  Pressure  Distribution 
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this  portion  of  a  sphere  is  not  the  same  as  that  over  a  sphere  alone  in 
incompressible  flow.  The  integration  of  the  boundary  layer  equations 
continues  from  the  stagnation  point  along  a  meridian  until  the  sphere- 
afterbody  interface  has  been  reached.  This  stopping  point  is  designated 
by  the  angle  attaining  a  particular  maximum  value.  This  maximum  value 
is  a  function  of  the  circumferential  position  on  the  interface  and 
can  be  calculated  from  Equation  (4.2).  Beyond  the  interface,  the 
inviscid  surface  streamlines  must  be  integrated  numerically. 
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SECTION  6 


DESCRIPTION  OF  COMPUTATIONAL  METHOD 

The  method  presented  traces  inviscid  surface  streamlines  while 
simultaneously  computing  the  properties  of  the  boundary  layer  up  to 
the  separation  point.  Tracing  a  streamline  involves  the  numerical 
integration  in  a  step-by-step  fashion  of  Equations  (2.9),  (2.10),  and 
(2.18)  to  determine  the  axial  position,  circumferential  angle,  and 
the  streamline  angle  (see  INVISCID  SURFACE  STREAMLINES).  In  conjunc¬ 
tion  with  the  differential  equations  for  the  streamlines.  Equations 
(2.31)  and  (2.37)  are  also  integrated  to  give  the  scale  factor  along 
the  streamlines.  After  each  integration  increment  along  a  stream¬ 
line,  the  boundary- layer  equations  are  then  integrated  by  either 
Hall's  or  Blotter's  method  to  determine  the  local  velocity  Drofile 
across  the  boundary  layer  (see  BOUNDARY  LAYER  METHODS).  This  pro¬ 
file  is  used  to  determine  the  local  value  of  Cf  KA  .  The  seoara- 
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tion  point  for  the  flow  along  a  particular  streamline  is  assumed  to 
occur  when  this  parameter  passes  through  zero.  All  calculations  stop 
at  this  point  since  both  the  streamline  and  boundary-layer  equations 
are  invalid  in  the  separated  region.  Several  streamlines  are  calculated 
to  get  a  distribution  around  the  body. 

In  order  to  begin  the  boundary  layer  integration,  it  is  first 
necessary  to  establish  the  initial  boundary  layer  velocity  profile  at 
the  stagnation  point.  For  Blottner's  method  this  neccesitates  solving 
the  similar  F-V  equations  while  for  Hall's  method,  each  point  in  the 
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profile  is  identically  zero.  For  instances  in  which  the  pressure  distri¬ 
bution  is  expressed  analytically,  the  integration  of  the  boundary  layer 
continues  along  an  inviscid  surface  streamline  in  increments  of  As 
(which  has  units  of  the  input  geometry).  For  cases  in  which  only 
experimental  pressures  are  available,  the  angle  \ p  is  first  calculated 
given  a  position  on  the  sphere-afterbody  interface.  The  integration  of 
the  boundary  layer  then  proceeds  in  a  step-by-step  fashion  in  increments 
of  the  angle  \p,  Ai|>,  on  the  spherical  cap  until  the  value  of 
1 1>  at  the  interface  has  been  reached.  The  boundary  condition  on  the 
fluid  velocity  at  the  edge  of  the  boundary  layer  is  a  function  of  ip  only 
and  is  obtained  at  each  step  during  the  Integration  from  subroutine 
SPHCAP.  The  integration  of  the  streamlines  begins  at  the  sphere- 
afterbody  interface. 

The  boundary-layer  profile  convergence  test  is  then  applied  after 
each  iteration  of  the  solution.  This  computer  program  employs  two 
options  with  which  to  define  a  converged  solution.  One  option  requires 
that  the  skin  friction  parameter,  |cf-/ReL.jw»  change  by  less  than  0.5 
percent  between  successive  iterations  in  order  for  the  solution  to  be 
considered  to  have  converged.  The  other  option  requires  that  the 
velocity  at  each  grid  location  change  by  less  than  0.1  percent  from 
the  previous  iteration. 

Once  the  solution  has  converged,  the  edge  test  is  oerformed.  This 
test  effectively  determines  whether  the  point  at  the  edge  of  the  boundary 
layer  spans  the  total  thickness.  There  are  two  options  regarding  the 
edge  test.  One  option  requires  that  the  velocity  at  the  grid  location 
just  inside  the  imposed  boundary -layer  edge  be  at  least  99.95  percent 


of  the  velocity  at  the  boundary  in  order  for  the  imposed  thickness  to  be 
adequate.  The  second  option  requires  that  the  skin  friction  parameter, 
Cf  ,at  the  edge  be  less  than  or  equal  to  0.005.  If  the  option 
employed  should  fail,  an  additional  point  is  added  at  the  boundary -layer 
edge  of  both  the  present  and  previous  computational  station.  The 
velocity  at  the  outer  grid  point  of  the  previous  station  is  assigned 
the  value  of  the  edge  velocity  at  that  station.  With  the  addition  of 
the  point  at  the  edge,  the  calculations  for  the  present  station  are  re¬ 
peated.  This  procedure  is  followed  until  both  tests  have  been  satisfied. 
If  the  number  of  points  added  at  the  edge  should  eventually  exceed  50, 
the  step  size  in  the  normal  direction  is  doubled  and  every  other  point 
within  the  boundary  layer  is  discarded.  At  this  time  the  step  size 
along  the  streamline  is  also  doubled. 

Three  methods  are  available  with  which  to  integrate  the  inviscid 
streamline  differential  equations.  The  predictor-corrector  method  of 
Milnes  (Ref.  12)  features  rapid  execution  and  has  been  incorporated  into 
the  computational  code.  This  method,  however,  is  not  self-starting 
and  makes  no  check  for  truncation  error  (see  SUBROUTINE  MILNES).  The 
method  of  Gear  is  used  to  generate  the  starting  values.  Gear's  method 
(Ref.  12)  is  useful  for  instances  in  which  a  stiff  system  of  first-order 
differential  equations  is  being  Integrated  (see  SUBROUTINE  DGEAR).  The 
last  method  is  the  fourth-order  Runge-Kutta  method.  This  method  is  not 
as  well  suited  to  stiff  systems  because  the  step  size  becomes  prohibi¬ 
tively  small  In  the  attempt  to  minimize  the  truncation  error  during 
Integration  (see  FUNCTION  KRUNGE). 


The  boundary  layer  is  integrated  in  increments  of  As  along  a 
streamline.  .The  boundary  condition  on  the  fluid  velocity  at  the  outer 
edge  is  obtained  from  subroutine  INVISD.  If  experimental  pressures  are 
supplied,  a  second  subroutine,  SPHCAP,  provides  the  necessary  condi¬ 
tions  for  points  on  the  spherical  cap.  The  integration  of  the  boundary 
layer  continues  up  to  the  point  at  which  fCf  jRg  |  reaches  or  passes 
through  zero. 

At  larger  angles-of-attack,  the  streamlines  auite  frequently  wrap 
around  the  body  so  rapidly  that  it  is  difficult  to  resolve  the  boundary 
layer  at  points  further  down  the  body.  A  technique  which  employs  a  shift 
from  the  windward  streamline  may  be  implemented  in  order  to  accomplish 
this.  The  integration  of  the  boundary  layer  continues  along  the  wind¬ 
ward  streamline  to  the  input  axial  position,  XMAX.  At  this  point  the 
circumferential  angle  is  changed  from  zero  to  one  degree.  From  this 
point  on  the  integration  continues  along  this  newly  defined  streamline 
to  the  separation  point.  With  this  technique  it  is  possible  to  trace 
streamlines  that  otherwise  would  have  been  unobtainable. 

The  separation  point  along  a  streamline  is  approximated  by  linear 
interpolation  using  the  last  two  converged  solutions  since  the  boundary- 
layer  Drofile  frequently  falls  to  converge  once  in  the  separated  region. 
This  entire  procedure  is  repeated  for  each  of  the  streamlines.  The 
total  number  of  streamlines  is  an  input  parameter  called  KBM. 
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SECTION  7 

RESULTS  AND  DISCUSSION 


In  order  to  illustrate  the  validity  of  the  techniques  employed  in 
the  computer  program,  results  are  presented  for  a  sphere,  ellipsoid  of 
revolution  at  an  angle  of  attack  and  a  sphere-ogive-cylinder  configura¬ 
tion  at  an  angle  of  attack.  Each  test  case  represents  a  step  up  in  the 
complexity  of  the  analysis.  The  results  for  each  geometry  consist  pri¬ 
marily  of  a  comparison  between  the  solutions  obtained  by  both  Hall's  and 
Blottner's  methods.  Additional  results  generated  for  the  ellipsoid  of 
revolution  at  two  different  angles  of  attack  are  compared  to  fully  three- 
dimensional  boundary- layer  calculations.  All  computations  were  performed 
on  the  IBM  370/165  digital  computer  at  North  Carolina  State  University. 
Computer  times  in  this  section  are  in  CPU  seconds.  All  cases  start  with 
40  points  across  the  boundary  layer. 


7.1.  Sphere 

The  sphere  geometry  provides  the  opportunity  to  validate  the  com¬ 
putational  code  itself.  The  comparisons  presented  in  Tables  3.1,  3.2, 
and  7.1  serve  as  verification.  Results  generaged  on  both  a  cylinder 
and  a  flat  plate  compared  quite  well  with  the  accepted  values. 

7.2.  Ellipsoid  of  Revolution 

The  ellipsoid  selected  for  this  case  had  a  thickness  ratio  of  1/4, 
a  total  length  L  =  2a,  and  was  examined  at  both  12°  and  30°  angle-of- 
attack  (see  Figure  7.1).  The  potential  solution  (Ref.  13)  was  available 
in  the  form  of  an  analytical  expression  (see  Appendix  A. 6).  For  this 
case,  only  the  differential  equations  in  (2.9),  (2.10)  and  [2  31)  were 
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Table  7.1.  Effect  of  Convergence  and  Edge  Criteria  on  Skin  Friction 
for  Sphere  in  Incompressible  Flow 


As  =  0. 

Hall 
05  An 

wmmm 

(jj  Blottner 

As  =  0.05  An  = 

0.11539 

<l» 

NC  =  0 

NC  =  1 

NC  =  0 

NC  =  1 

20. 

2.578 

2.579 

2.710 

2.573 

40. 

2.430 

2.418 

2.431 

2.417 

60 

2.150 

2.149 

2.150 

2.149 

80 

1.656 

1.656 

1.654 

1.657 

90 

1.249 

1.276 

1.271 

1.271 

100 

0.650 

0.645 

0.640 

0.642 

104 

0.202 

0.193 

0.226 

0.226 

integrated  since  analytical  expressions  for  the  streamline  angle.  Equa¬ 
tion  (2.18),  and  its  circumferential  derivative.  Equation  (2.39),  were 
available  (see  Appendix  A. 7).  For  these  cases.  As  =  0.05  and  An  =  0.0471 
were  used  for  Hall's  method,  and  As  =  0.05  and  An  =  0.115385  for  Blottner's 
method. 

The  results  of  both  Hall's  and  Blottner's  methods  are  presented  for 
a  variety  of  streamlines  on  this  configuration  in  Table  7.2.  This 
table  includes  results  that  were  obtained  from  the  streamline  shifting 
technique  described  in  the  section  of  this  thesis  labeled  DESCRIPTION 
OF  COMPUTATIONAL  METHOD.  XMAX  is  the  axial  position  at  which  a  stream¬ 
line  shift  from  a  circumferential  position  of  zero  to  one  degree  was 
made.  The  results  of  both  methods  agree  quite  well  despite  the  fact 
that  the  step  size  along  a  streamline  differs  between  the  two  methods 
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Table  7.2.  Comparison  of  Separation  Points  Between  Hall's  and 
Blottner's  Methods  on  Ellipsoid  of  Revolution  with 
Thickness  Ratio  1/4  at  30°  Angle-of -Attack 


Method 

Beta 

XMAX/L 

Points 

Added 

Computer  Separation 

Time  (sec)  X/L  <j> 

Hall 

20.0 

-- 

52 

6 

0.209 

143.28 

50.0 

-- 

55 

7 

0.222 

142.76 

0.0 

0.30 

53 

13 

1.363 

104.42 

0.0 

0.50 

49 

13 

1.626 

95.89 

0.0 

0.80 

48 

14 

1.857 

82.18 

Blottner 

20.0 

— 

3 

7 

0.212 

143.94 

50.0 

-- 

3 

7 

0.220 

142.43 

0.0 

0.30 

1 

12 

1.366 

104.15 

0.0 

0.50 

2 

12 

1.640 

97.46 

0.0 

0.80 

0 

12 

1.857 

82.16 

by  a  factor  of  20.  Hall's  method  generally  required  a  greater  amount 
of  computational  time.  This  is  most  likely  due  to  the  greater  number  of 
points  that  had  to  be  added  at  the  imposed  boundary  layer  edge. 

Tables  7.3  and  7.4  compare  results  generated  by  the  axisymmetric 
analogue  using  Hall's  and  Blottner's  methods  to  the  three-dimensional 
boundary  layer  calculations  of  Wang  (Ref.  13).  In  both  cases  the  com¬ 
parison  is  increasingly  degraded  as  the  leeside  of  the  body  is  approached. 
In  this  instance,  both  tables  suggest  that  the  axi symmetric  analogue 
yields  quite  good  results  on  the  windside  of  the  body.  The  separated 
region  for  this  case  is  shown  graphically  in  Figure  7.2.  The  separated 
region  as  calculated  by  Wang  is  also  included  for  comparison. 
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Table  7.3.  Comparison  of  Separation  Points  Between  Hall's  Method 
and  Three-Dimensional  Boundary  Layer  Calculations  on 
Ellipsoid  of  Revolution  with  Thickness  Ratio  1/4  at  3f)° 
Angle-of-Attack 


Circumferential  Separation, 

<|>,  (Degrees) 

Axial  Station  (X/L) 

Hall 's  Method  3- 

-D  Results 

0.209 

143.28 

131.25 

0.222 

142.76 

130.00 

0.328 

132.56 

125.50 

1.363 

104.42 

102.50 

1.626 

95.89 

95.00 

1.857 

82.18 

82.50 

Table  7.4.  Comparison  of  Separation  Points  Between  Blottner's  Method 
and  Three-Dimensional  Boundary  Layer  Calculations  on 
Ellipsoid  of  Revolution  with  Thickness  Ratio  1/4  at  30° 
Angle-of-Attack 


Circumferential  Separation,  4>,  (Degrees) 

Axial  Station  (X/L) 

Blottner's  Method  3-D  Results 

0.212 

143.94 

131.25 

0.220 

142.43 

130.00 

1.366 

104.15 

102.50 

1.640 

97.46 

94.50 

1.857 

82.16 

82.50 

110° 
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Figure  7.2.  Separation  Region  on  Ellipsoid  of  Revolution  with  Thickness 

Ratio  1/4  at  30  Degrees  Angle-of-Attack 


Figure  7.3  depicts  the  variation  of 


Re^  in  the  windward  plane 


'  'OJ 

for  the  axisymmetric  analogue  using  Hall's  method  and  Wang's  fully  three- 
dimensional  approach  for  an  angle-of-attack  of  12°  (Ref.  5).  The  two 
methods  compare  reasonably  well. 


7.3.  Sphere-Ogi ve-Cyl inder 

This  configuration,  whose  geometry  is  depicted  in  Figure  5.1,  was  in¬ 
put  to  the  program  in  dimensionless  form.  The  normalizing  quantity  was 
the  cylinder  radius  which  measured  3.8  inches.  The  total  non-dimensional 
body  length  was  17.5.  This  geometry  was  investigated  at  45°  angle-of- 
attack.  The  experimental  pressure  data  consisted  of  discrete  pressure 
coefficients  distributed  along  30  axial  stations,  each  having  10  circum¬ 
ferential  stations.  Pressure  data  for  the  region  0  _  x  13.75  was  cal¬ 
culated  by  the  USSAERO  panel  method  while  the  pressure  over  the  section 
0.92  1  x  ^  17.5  consisted  of  actual  experimental  pressures  obtained  from 
the  wind  tunnel.  Before  implementing  these  pressures  into  the  computa¬ 
tional  code,  it  was  necessary  to  smooth  and  interpolate  the  data  in  the 
region  where  the  pressures  overlapped.  Interpolated  data  were  used  to 
form  additional  axial  stations  near  the  nose  since  large  pressure  grad¬ 
ients  are  present  on  the  forward  portion  of  the  body.  This  required 
the  addition  of  four  more  axial  stations  in  that  region  (the  resulting 
pressure  coefficients  as  well  as  the  remaining  program  inputs  are  pre¬ 
sented  in  Appendix  A. 22). 

The  technique  developed  to  represent  the  pressure  distribution  on 
the  spherical  cap  by  a  10- term  Fourier  cosine  series  was  found  to  perform 
only  marginally.  The  series  Drovided  continuity  in  the  pressure  coefficient 
across  the  interface  (when  the  quadratic  spline  was  first  employed)  but 
did  not  necessarily  provide  continuity  in  the  related  derivatives.  In  instances 
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Wang's  (5) 


riction  Distribution  in  Windward  Plane  of  ElliDsoid  of  Revolution 
Thickness  Ratio  1/4  at  12  Degrees  Angle-of-Attack 


In  which  the  pressure  derivatives  were  not  continuous  across  the  inter¬ 
face,  the  one-dimensional  quadratic  spline  technique  was  extended  from 
the  interface  to  the  stagnation  point  to  model  the  pressure  variation 
on  the  spherical  cap  (see  SURFACE  PRESSURE  DISTRIBUTION).  The  resulting 
pressure  variation  was  a  function  of  the  circumferential  position  on 
the  interface.  Despite  this,  the  streamlines  in  this  region  were  still 
assumed  to  follow  spherical  meridians. 

The  doubly  quadratic  spline  technique  employed  in  this  computer 
program  was  found  to  model  the  pressure  coefficient  variation  quite 
satisfactorily.  This  technique  requires  that  the  axial  derivative  of 
the  pressure  coefficient  at  the  interface  and  body  end  for  each  cir¬ 
cumferential  plane  be  known.  While  the  pressure  coefficient  across  the 
interface  was  continuous,  in  most  cases  the  axial  derivative  was  not 
and,  hence,  was  also  supplied  as  program  input  (rather  than  calculated 
in  the  program). 

The  streamline  angle  and  its  circumferential  derivative  were  cal¬ 
culated  in  this  case  by  numerically  integrating  Equations  (2.18)  and 
(2.39).  These  equations  are  functions  of  the  inviscid  edge  velocity, 
the  pressure  coefficient  and  its  derivatives.  These  parameters  were 
provided  by  the  spline  fit.  Although  the  calculated  second  derivatives 
of  the  pressure  coefficient  in  the  circumferential  direction  are  con¬ 
stant  in  the  interval  in  which  the  quadratic  is  used,  they  were  found 
to  be  accurate  enough  to  be  used  in  the  integration  of  Equation  (2.39). 
Table  7.5  provides  a  comparison  between  the  calculated  separation 
points  using  Hall's  and  Blottner's  methods.  The  corresponding  separa¬ 
tion  points  agree  very  well.  Note  that  for  the  cases  of  B  =  145"  and 
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Table  7.5.  Comparison  of  Separation  Points  Between  Hall's  and 
Blottner's  Methods  on  Sphere-Ogive-Cylinder  Config¬ 
uration  at  45°  Angle-of-Attack 


Method 

Beta 

Points 

Added 

Computer 
Time  (sec) 

X/Rc 

Separation 

<P 

Hall 

50 

39 

9 

0.2158 

119.15 

-- 

75 

0 

8 

0.1777 

125.32 

-- 

100 

3 

8 

0.1562 

121.34 

-- 

145 

5 

8 

-- 

-- 

115.14 

160 

27 

8 

-- 

-- 

119.64 

Blottner 

50 

3 

8 

0.2159 

119.21 

-- 

75 

3 

8 

0.1763 

123.75 

-- 

100 

1 

8 

0.1564 

121.66 

-- 

145 

2 

8 

-- 

-- 

115.16 

160 

0 

8 

-- 

-- 

119.25 

160°,  the  flow  separates  while  on  the  spherical  cap.  If  the  pressure 
distribution  had  been  truly  axisymmetric,  the  angle  of  separation,  \p, 
would  have  been  identical  in  each  case.  The  computational  time  required 
for  this  configuration  was  greater  than  that  for  the  ellipsoid  of 
revolution  though  still  quite  reasonable. 

Information  relative  to  the  step  sizes  and  spacings  is  given  in 
Table  7.6.  The  step  size  on  the  spherical  cap  was  Aip  =  2°  in  each  case. 
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Table  7.6.  Computational  Step  Sizes  and  Spacings  for  Results  in 
Table  7.5 


Method 

Beta 

As 

An 

or  An 

Steps  to 
Separation 

Hall 

50 

0.001 

0.003846 

109 

75 

0.004 

0.015385 

50 

100 

0.004 

0.007692 

52 

145 

— 

0.007692 

160 

— 

0.003846 

55 

Blottner 

50 

0.001 

0.1153846 

136 

75 

100 

145 

— 

0.1153846 

59 

160 

— 

0.1153846 

60 
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SECTION  8 


CONCLUDING  REMARKS 

A  method  is  developed  for  calculating  laminar  boundary  layers  along 
inviscid  surface  streamlines  on  axisymmetric  bodies  at  angles  of  attack 
in  incompressible  flow.  By  application  of  the  axisymmetric  analogue 
concept  in  the  present  technique,  a  substantial  savings  in  computer  time 
over  fully  three-dimensional  boundary  layer  techniques  may  be  realized. 

The  boundary  layer  integration  techniques  of  Hall  and  Blottner  were 
found  to  compare  exceptionally  well  with  each  other  on  each  of  the 
geometries  investigated.  Results  generated  on  the  windward  plane  of  an 
ellipsoid  of  revolution  with  thickness  ratio  1/4  and  angle  of  attack  of 
12°  compared  satisfactorily  with  results  generated  by  a  fully  three- 
dimensional  technique.  The  separation  points  calculated  by  the  present 
technique  for  a  variety  of  streamlines  on  the  same  ellipsoid  of  revolu¬ 
tion  at  30°  angle  of  attack  were  in  fair  agreement  to  those  generated 
by  a  three-dimensional  technique.  The  comparison  was  generally  better 
on  the  windside  of  the  body. 

The  series  expression  used  to  model  the  pressure  coefficient  on  the 
spherical  cap  of  the  sphere-ogi ve-cylinder  configuration  performed  only 
satisfactory.  The  technique  preserved  continuity  in  the  pressure  co¬ 
efficient  across  the  interface  (the  point  at  which  the  quadratic  spline 
technique  was  implemented)  but  did  not  provide  continuity  of  the  axial 
derivative.  To  circumvent  this  problem,  the  one-dimensional  quadratic 
spline  was  extended  to  the  stagnation  point  on  the  spherical  cap.  De¬ 
spite  this,  the  assumption  of  spherical  streamlines  was  still  made  with 
reasonable  accuracy.  The  doubly  quadratic  spline  renresentation  of  the 
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pressure  coefficient  on  the  afterbody  was  found  to  perform  quite  well  as 
long  as  an  adequate  number  of  pressure  stations  were  input.  Oscillations 
in  the  pressure  function  were  generally  less  frequent  than  might  be 
expected  if  other  techniques  had  been  used. 

The  relative  inexpense,  coupled  with  reasonable  accuracy  makes  the 
present  method  attractive  for  preliminary  design  studies.  Further  com¬ 
parisons  with  fully  three-dimensional  boundary  layer  calculations  are 
necessary  in  order  to  more  thoroughly  evaluate  the  apDl icabi 1 i ty  of  the 
axisymmetric  analogue  in  subsonic  flow. 
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LIST  OF  SYMBOLS 


coefficients  defined  in  Equation  (4.3) 

coefficients  for  finite-difference  boundary-layer 
equations,  defined  in  Equations  (3.8),  (3.9),  (3.12), 
(3.13) 

unit  vectors  on  body  surface  along  body  meridian 
given  by  Equation  (2.2) 

unit  vectors  in  streamline  coordinate  system  given 
by  Equations  (2.3),  (2.4)  and  (2.1) 

unit  vectors  in  cylindrical  coordinate  system,  (see 
Figure  2.1) 

2t 

skin  friction  coefficient,  — “ 

(K 

total  derivative  along  streamline 

ratio  of  local  velocity  to  velocity  at  boundary-layer 
edge,  as  defined  by  Equation  (3.6),  dimensionless 

scale  factor  in  3  direction,  dimensionless 

arbitrary  constant  for  Equation  (3.1)  and  (3.2) 

body  length,  dimensionless 

coordinate  normal  to  body  surface  and  streamline 
dimensional  pressure,  lb/ ft7  or  N/m2 
body  radius,  dimensionless 

radius  of  cylinder  in  sphere-ogive-cylinder  con¬ 
figuration,  dimensionless 

radius  of  spherical  cap,  dimensionless 

puj. 

freestream  Reynold's  number, 

distance  along  streamline,  dimensionless 

local  fluid  velocity  in  boundary  layer  (in  direction 
of  a  streamline)  ft/sec  or  m/sec 


u 

v 

Y 

x,y,z 

Y(l) 

Y(2) 

Y  ( 3 ) 

Y  ( 4 ) 

Y  ( 5 ) 

Y  ( 6 ) 

a 

8 

6 

r 

5 

n 

4> 

6 

<P 

p 

v 


inviscid  fluid  velocity,  ft/sec  or  m/sec 

local  fluid  velocity  normal  to  streamline  and  body 
surface ,  ft/sec  or  m/sec 

parameter  defined  by  Equation  (3.5) 

body  geometry  coordinate  axes  (see  Figure  2.1) 

axial  position,  x,  dimensionless,  ft  or  m 

circumferential  angle,  d> ,  rads 

streamline  angle,  0,  rads 

90 

34>  x 
In 

aex 

transformed  streamline  coordinate,  r 
angle  of  attack,  degrees 

coordinate  normal  to  streamline  and  tanqent  to  body 

2F  due 

pressure  gradient  parameter,  ~ 

body  angle,  radians  (see  Figure  2.2) 

transformed  streamline  coordinate  as  defined  in 
(3.1) 

transformed  coordinate  normal  to  body  surface 
defined  in  (3.2) 

circumferential  angle  (see  Figure  2.1),  rads 

streamline  angle  (see  Figure  2.3),  rads 

angle  between  stagnation  line  and  radius  vector  (see 
Figure  4.1),  rads 

density,  slug/ft3  or  kg/m1 

coefficient  of  viscosity,  sluq/ft-sec  or  kg/m-s 


64 


Subscripts 

e 

eff 

0) 

0 

i 

I 

j 

SP 


frees tream  conditions 

edge  of  the  boundary  layer 

effective  value 

at  the  wall 

initial  value 

streamline  grid  index 

value  at  sphere-afterbody  interface 

normal  grid  index 

value  at  stagnation  point 


Superscripts 

_  denotes  a  quantity  from  a  previous  iteration 

*  denotes  a  dimensionless  quantity 
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APPENDIX  A.  EQUATIONS 
1«  Derivation  of  Fquations  (3.3)  _and_ J 3.4) 

Blottner's  method  involves  the  boundary  layer  equations  written  in 
F-V  similarity  form.  These  are  obtained  by  application  of  the  Lew- 
Lees  transformation  which  is  defined  as 


and 


£(s)  =  Kppu^ 


ds 


(A. 1.1) 


u  prn 

n{s,n)  =  — -  fK  (A. 1.2) 

/W. 


for  incompressible  flow.  The  transformation  operators  may  be  con¬ 
structed  and  expressed  as 


9  _  ^e  2  _8_  t  8n  9 

3s  U  9f  9s  9r) 

no  ' 


(A. 1.3) 


and 


a  _  uepr  1  8 

3n  "  m  /pph7  9ri 


(A. 1.4) 


where  the  arbitrary  constant,  K,  has  been  assiqned  the  value 


K  =  1/ppu 

m 


(A. 1.5) 
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The  dependent  variable  in  Blottner's  equations  is  defined  as 


F  -  u/ue  .  (A. 1.6) 

Application  of  each  transformation  operator  and  the  definition  of  F  to 
the  continuity  equation. 


i(£ui  + 

3s 


(A. 1.7) 


yields 


^rueF>  „  an 
3 f  ~  3s  rue 


IE 

3ri 


V^p 


(A. 1.8) 


Note  that  neither  uQ  nor  r  are  functions  of  the  normal  coordinate  n. 
Expanding  this  equation  yields 


—  r  \ru 


EE 

3f 


c  dr 

UeF  -rr  + 


dr 


fr 


du 
_ e 

dr. 


a  3n  3F  . 
+  ru.  ~  -r-  + 


u  r2p 
e 


e  85  3ri  m  /pi^r 


3 v 
3n 


TT  =  0 


(A. 1.9) 


This  equation  may  be  rewritten  as 


"ki  k  ♦  f  dr  t  .L  +  ru  5a  if.  t  vk 

u,„  j3^  r  &'■  ug  dr  j  e  3s  3n 


pyu 


Is  -  0 

3n 


(A. 1.10) 


By  application  of  the  product  rule,  this  equation  may  be  written  as 


3F  +  F  dr  +  X  +  .  U<"  JL 

3f  r  dr,  u„  dr  uTr7  3n 
e  e 


u  r'pv 


|n  F  + 
e  8s  /2T  /Jm 


_U,X’  r  k  3ri) 

3l)  gs 


J 


(A. 1.11) 
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The  derivati/e  in  the  last  term  of  this  equation  may  be  rewritten  as: 


_3_ 

[3n] 

_  _3_  3n 

f3nl 

3n 

3s; 

3n  3n 

[as 

/Jt,  / PMUr> 


rn _ 

/  2f,  J  ppiT 


+ 


dr 

ds 


2  r. 


uerp 


/jrr 


(A. 1.12) 


which  becomes,  upon  making  use  of  Equations  (A. 1.1)  and  (A. 1.5), 


JL  fin' 
3n  ,3s 


2  du  u  r  .  u  r2 
r_ _ e  _e_  dr  _  e 

u  d r,  u  dr,  '  2u  r 


Substitution  of  this  relation  into  Equation  (A. 1.11)  yields 


3F  F  dr  F  due 
3C  r  df,  ue  dC 


3  jy*1  in  c  + 
3n  liTr7  3s 


u„pv 

uer  m 


JL^e  _  Fdr 
ug  dc  r  d £ 


(A,  1.13) 


This  equation  may  then  be  written  as 


25  f  *  F  + 


r=  0 

3n 


(A. 1.14) 


where 


V  =  2  £ 


u 


3n 
3s  \TP 


PVKn _ 1 

/  2£  / ppu  u  r 
,r>  e 


(A. 1.15) 


Application  of  the  transformation  operators  in  (A. 1.3)  and  (A. 1.4) 
and  the  definition  of  F  to  the  momentum  equation 


(A. 1.16) 


du 


3u  ,  3u  _ 

u  3s  v  3n  ~  ue  ds  '  p  3n 


e  +  li  32u 


yields 


“eF 


f.y-'.Ku/)  m 

u  3r  3s  e  3n  | 


vu2rp 

e 


.  if. 
/If.  /mT  3n 


u2r2  du 
_e _ e 

u  dr 


uer  p  i-  32F 
2r,puu  f>  Sn'2 


(A. 1.17) 


This  may  be  expanded  to  give 


(u  rF)?  du  u  'r2r  ~c  ~  »c 

_e _ _ _ e  _e _  3F  2f  3n  3F 

U  dr,  u  3r  V  3s  3n 


vue  rp  3F 
3n 


/If,  /pvuT 


(u  r)2  du 
e  _ e 

u_  dr 


uer?  3?F 
2ru  3 n7 


(A. 1.18) 


This  equation  may  then  be  rearranged  to  yield 


0,  du  ,.r  2f,Fu  _  »r  w  2f,u 

2f.  r?  e  ,  9fr  3F  _ 3r|  3F  v _ « 

u  dT  '  ar  u  r7  3s  3q  rrr  u_r 
e  e 


_£ _  3_F 

/IT  uer  JmT  3n 


2T  ^e  32F 
Ug  dr,  Srp" 


(A  .1.19) 


or,  finally. 


3F 


3F  3?F 


2TF  ~  +  (F2  -  1)  ft  +  V  £  -  24-  =  0 


ar 


3n  3ri2 


(A. 1.20) 


where  p,  is  the  pressure  gradient  parameter  and  is  defined  as 


or  du 
2f  _  e 

Up  dr' 
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2.  Derivation  of  Equations  (3 .8)  and  (3.9) 

At  all  points  off  the  stagnation  point,  the  full  system  of  F-V 
equations  must  be  solved.  The  continuity  equation. 


2F  +  |^  +  F  =  ° 

3F.  3n 


(A. 2.1) 


is  evaluated  at  the  point  (i  +  1/2,  j  -  1/2)  and  may  be  expressed  as 


3F 

3C 


3F 

3C 


i+H.M  .  3 V 
3n 


i+’sj-'i 


4 


=  0 


Substituting  second-order  accurate  expressions  for  the  appropriate 
quantities  yields 


( 

AF 


Fi+l,j  ”  Fi , j  +  Fi+l,j-l  "  Fi,j-1 


+  ^Vi+I2».i  '  Vi+12,j-l) 
An 


+  +  *  Fi.o  _+  F:UL-ll  =  0 

4  u 


After  rearranging,  the  continuity  equation  may  be  expressed  as 


A2Fi  +  l,j-l  +  B2Fi  +  l,j  +  C2Vi+'„j-l  +  E2Vi+'.,.i  =  °2  ^.2.2) 
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where 


A2  =  An('4  +  AfJ 
B2  =  An(‘4  +  Ci+1?/Ar) 


and 


The  momentum  equation, 

2fF  I*  V  -  »  -  0“  0  (A-2-3) 

is  evaluated  at  (i+'2,j).  The  terns  of  Equation  (A. 2. 3)  become,  resDectively, 


2rF  — - 
■  DC 


1+,a.j 


2;'i+,z 


Fi±l>j  *  Fi  ,j 


Fi  +  l.j  ~  FiJ 


AC 


DF 

=  vli 

DF 

+  i>r: 

Dr, 

i+l20 

2 

Dn 

i+1,. 

i  "n 

i  J 

*(F?  ~  ^i+CJ  =  » 


F’  .  +  F?  . 

o  A 


;>?F 

a»r 


i+CJ 


D?F 

Dr,7 


+  3fF| 

3ri  lij 


(A. 2. 4) 
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Linearizing  the  first  three  terms  using  the  Newton- Raohson  method  yields 

2f  F  —  =  F  -  F2  _  c2  ) 

i+4  i+^.j  Af,  Uhi+l,j  i+l,j  hi+l,j  1  ,j  ’ 


v  ?1 

3n 


=  —  V 

i+w  2  1n-J 


3F 

+  It 

3n 

i  +  1 ,3 

3n 

i  ,J 

+  1 V  f3-F-l 

2  V',,j  [3nJ 


h+1,3 


1  v  i£| 

2  i+'s, j  an 


i+lj 


®<F’  -  ■  ~z' 


a  (2F.a1  .  F.  .  .  -  F2  .  +  F?  .  -  1)  .  (A. 2. 5) 

1+1,3  1+1.3  1+1.3  1,3 


Substituting  second-order  accurate  finite  difference  expressions  in  the 

second  of  Equations  (A. 2. 5)  and  the  last  of  Equations  (A. 2. 4)  gives, 
respectively. 


3F| 


r,  -  f,  , 


+  —  v 

2  i+%,3 


i+'a.j 

An 

An 

Fi+1 ,3+1  ’  Fi+1,3-1 

Fj  +  l,j+_l  '  F i + 1  ,j - 1 

2An 

2  i+',,3 

2Ar, 

J 

and 


32F I 


i+?2,j 


Fi+i,j+l  ~  2Ft+l,J. 


A>r 


Fi  +  l,j-l  .  F i_, j +J_ 


2F .  . 
—bJ. 
Air 


+  F. 


After  substituting  the  appropriate  expressions  into  Equation  (A. 2. 3), 


solving  for  the  unbarred  quantities  and  rearranging  ,the  equation  may  be 
written  as 


AlFi  +  l,j-l  +  BlFi  +  l  ,j  +  ClFi  +  l,,i  +  l  +  ElVi  +  ’,.,i  =  D1 
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where 


A1  ■  -  ?  <>  +  2-  ^  V„j> 

B,  =  1  ♦  An'  Ft+lo.  (F)t,  * 


C1  "  2  ^  '  2  An 


E.  =  4-  An  (F.  ...  -  F.  .  .  +  F. . . 

1  4  x  l  ,j+l  i  ,j-l  i  +  l» 


.i+l 


and 


D1  ’  I  (Fi,.1+l  '  2Fi,j  *  *  1  An2  Sit'. 


(1  +  F? . .  . 
i+l, J 


Ml  -  F-  .) 

l,.l 


+  4  Ar|  Vi+',,.i  (Fi  +  l,j  +  l  "  Fi  +  l,j-l) 


*  *  V*  <Fm.j  * ri,.i)/Af' 
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3.  Derivation  of  Fini te-Difference  F-V  Similarity 
Equations  at  Stagnation  Poin t 

At  the  stagnation  point  r  =  0  and  the  F-V  equations  reduce  to 

In  +  F  =  0  (A  -3.1) 

for  the  continuity  equation  and 

+iT<r2  -  ‘I  -0-0  (A -3.2) 

for  the  momentum  equation. 

The  continuity  equation  is  evaluated  at  (.i  -  1/2)  and  may  be 
exDressed  as 


(Fi  +  F 


=  0 


This  may  be  rewritten  as 


vj  ■ 


j-1 


Ari 

2 


(Fi  *  fm> 


(A  .3.3) 


The  momentum  equation  is  evaluated  at  (j)  after  first  being 
linearized  using  the  Newton-Ranhson  technique.  The  equation  mav  then 
be  written  as 


3F 

9f| 


+  V 


3_F 

3q 


-  V 


+  p  (2FF  -  P  -  1)  - 


;)?F 
.  ‘  '** 
:>rr 


=  0 
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where  the  barred  quantities  denote  the  expressions  from  the  previous 

iteration.  Substitutinq  appropriate  finite  difference  approximations 
for  each  of  the  above  terms  yields 


Solution  of  this  equation  for  the  unbarred  quantities  yields 


V.i-1  *  Vj *  V.m  +  Eiv.i '  Di  ■ 


where 


4.  Derivation  of  Equations  (3.10)  and  J3 . 1 1 J 

Hall's  method  involves  solving  the  boundary  layer  equations  written 
in  terms  of  dimensionless  primitive  variables.  The  transformations  used 
to  obtain  these  equations  are 

s*  =  s/L 

n*  =  f~RT  n/L 
CL 

u*  =  u/u 

or> 

v*  =  /TL-,  v/u 
L 

r*  =  r/L 


and 


and  a  star  denotes  a  dimensionless  quantity.  Accordingly,  the  trans¬ 
formation  operators  are  constructed  as  follows 


3  =  1  3_ 

3s  L  3s* 


and 
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Application  of  each  transformation  operator  to  the  continuity  equation, 


8s 


r 


8v 

8n 


=  0 


yields 


i'll  31X1^1  +  r*L  /~^L  _J_ 
L  8s*  L  8n* 


u  v” 


=  0 


which  will  simplify  to 


3( r*u*) 
8s* 


0 


Note  that  r  is  not  a  function  of  the  normal  coordinate  n. 

Application  of  each  of  the  transformation  operators  to  the  momentum 
equation. 


u 


Xu  +  v  = 

9s  8n 


yields 


8u* 

9s* 


vV 

nr* 


ARp.  -  .  u l  du  *  Rp.  * 

_“_L  9uX  _  _J*_’  *  e  p  _£l  8  u* 

L  9n*  L  ue  ds*  p  u<«  L  8n*z 


which  may  be  simplified  to 


u* 


du  * 
e 

ds* 


82u* 

8n*7 
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5.  Derivation  of  Equations  (3.12)  and  (3.13) 

The  nondimensional  continuity  equation, 

-^r1 +  "  £ ■ 0 

3s  3n 

(the  stars  have  been  deleted  for  clarity)  is  evaluated  at  the  point 
(i  +  1/2,  j  -  1/2)  in  the  computational  grid.  The  equation  may  then 
be  written  as 


1 

3(hu) 

,  3 ( hu ) 

+  h  ^ 

2 

3s 

i+52»j  35 

1+*Stj-l 

i+’z  3n 

=  0 


substituting  second-order  accurate  expressions  for  the  appropriate 
terms  yields 


1 

2 


(hu) . 


m,.i  ~  (hu)u  t  (hu)Hi.i-i  -  (hu)i.i-i 
As  As 


+  h 


i+,2 


An 


=  0 


Solving  for  v..,  .  gives 

I  +  2»  J 


.•  =  V 


An 


i+13,j  i+!j,j-l  2  Ash. 


1+,2 


h-,i  (u.,,  •  +  u.,.  .  .)  -  h.  (u.  .  +  u. 

i  +  l  i+l,.)  i  +  l,.)-r  i  i,.i  i , 


After  rearranging,  the  continuity  equation  may  be  rewritten  as 


Vi  +  lJ  *  B2ui  + 1 , j-1  +  C2vi+'„.)  +  r2vi+’,,.i-l  ~  D2 
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where 


2Ashj+s 

Anh1+1 


Anh . 

°2  =  2Ash7|7  (ui,J  +  Ui,M}  ‘ 


The  nondimensional  momentum  equation. 


9u  ,  9u  aue  ,  92u 
Os  v  9n  ue  ds  9n2 


is  evaluated  at  the  point  (i  +  1/2, j).  Taken  term  by  term,  this  may  be 


written  as 


u.,,  .  -  u.  . 

T  +  1..1  1.-1 


9u  _  j.  9u  9u 

*  tn.j  "  2  1*l,j  9ni,j 
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32u 

1 

32u 

♦  -A: 

3n2 

i+!2»j 

2 

3n^ 

eid 

i  ,j 

Substituting  second-order  accurate  expressions  for  the  remaining 
derivatives  yields 


3u 

3n 


1 


i+h  ,j 


4An 


Ui+l,j+l  '  ui  +  l,j-l  +  ui,j+l  ‘  ui,j-l 


and 


32u 

W 


i+1s>  j 


1 

2AnJ 


u.  ,  .  ,  -  2u.  ,  .  +  u.  ,  .  ,  +  u.  .  , 
i+l,  j+1  i+l, J  i+l, .1-1  i,j+l 


-  2u.  .  +  u.  .  , 
i,J  i,j-l 


Linearizing  using  the  Newton- Raphson  method  yields 


3u 

’  ’ 

u. ,  t  . 

1  +  1, .1 

3s 

i+ls,j 

As 

Ui+l,j  2As  ^ui  +  l,j  +  Ui  ,.i^ 


and 


3u 
v  3n 


=  Ity  /  77 


i+12,  j 


4An  ^Ui  +  l,j  +  l  Ui  +  l,.j-l  "*  Ui ,  j  +  1  U  i , .  j  - 1  ^ 
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After  solving  for  the  unbarred  quantities,  the  equation  may  be  written 
as 


Alui+l,j-l  +  BlUi  +  l  ,j  +  ClUi  +  l,.)+l  +  ElVi+1a. 
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6.  Derivation  of  Potential  Solution  for  Ellipsoid  of  Revolution 


In  cartesian  coordinates,  the  velocity  components  in  the  x,  y,  and 
z  directions  may  be  expressed  as  (Ref.  3) 


ux  = 


2u_  cos  a 


I?  2u  sin  a 
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(see  Figure  7. 5.  The  parameters  AQ  and  Bq  are  related  to  the  eccentricity 
and  are  defined  as  follows 


and 


1  +  e 


1  -  e 


-  e 


„  1  .  CL-  e!) 


2e 


£n 
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whe  re 

e  =  /(I  -  (b/a)2)  . 

In  cylindrical  coordinates,  the  velocity  components  may  be  expressed 
as 
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Substitution  of  the  appropriate  parameters  yields 
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The  component  of  the  total  velocity  along  a  body  meridian,  u  ,  may  be 
obtained  from 


um  .  -  u/ 


Substitution  of  the  appropriate  velocities  yields 


um 


2u  cos  a 

2u  sin  a 
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The  edge  velocity,  ug,  may  be  obtained  from 


u  =  /u  2  +  u . 2  . 

e  m  c|> 

Substitution  of  the  appropriate  expressions  for  the  velocity  components 

u  and  uA  then  gives 
m  cp 
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J.  Derivation  of  Analytical  Expression  for  Streamline  Angle 
and  Circumferential  Derivative  for  Ellipsoid 
of  Revolution 


The  streamline  angle  is  related  to  the  circumferential  and  meridianal 
velocity  components  and  may  be  shown  to  be 


Substitution  of  the  appropriate  velocity  components  (from  Appendix  A, 
Section  6)  yields 
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APPENDIX  B.  INPUT  PARAMETERS  AND  SUBROUTINES 
1.  Description  of  Input  Parameters 

Required  inputs  to  the  main  program  consist  of  the  following 
parameters  (note  that  all  parameters  describing  the  body  geometry  must 
be  input  with  the  same  units). 

RPER  radius  of  spherical  cap  (dimensionless,  ft  or  m) 

XNOSE  distance  from  body  nose  to  origin  of  body  axes 
(see  Figure  4. 1) (dimensionless ,  ft  or  m) 

BL  body  length,  L  (dimensionless,  ft  or  m) 

XINT  axial  location  of  sphere-afterbody  interface 
(see  Figure  4. 1) (dimensionless ,  ft  or  m) 

DST  maximum  step  size  along  a  streamline.  As  (dimensionless, 

ft  or  m) 

DPSI  step  size  in  degrees  of  arc  on  spherical  cap.  Ax  (degrees) 

ALPD  effective  angle-of-attack  ,  (degrees) 

NBS  Boundary  Layer  Method 

0  for  Hal  1 's  method 
1  for  Blottner's  method 

NT  Edge  Test 

0  for  point  added  at  edge  when  velocity  at  point  next 
to  edge  is  less  than  99.95X  of  point  at  edge 
1  for  point  added  when  (Cf /Rp^)^  is  above  D.005 
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NC  Convergence  Test 

0  for  convergence  based  on  (C^  /Re  changing  by 
less  than  0.57,  between  successive  iterations 
1  for  convergence  based  on  each  point  in  profile 
changing  by  less  than  0.1% 

MOT  Method  to  Integrate  Streamlines 

1  Milne's  predictor-corrector  method 

2  Runge-Kutta  method 

3  Gear's  method  for  stiff  system  of  differential 
equations 

MAXS  number  of  stations  to  be  computed 

KP  N  for  velocity  profiles  printed  every  Nth  station 

0  for  no  velocity  profiles  printed 

KPH  1  for  iterative  profiles  printed  after  each  Nth 

station 

0  for  no  iteration  profiles  printed 

KPO  Print  Out  Type 

0  for  ordinary  print  out 
I  for  additional  print  out 

KBM  number  of  streamlines  to  compute 

KBMS  Indicator  for  Streamline  Shifting 

0  for  ordinary  run  in  which  a  circumferential  position 
at  the  interface  is  specified 
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1  for  a  streamline  shift  to  be  made  from  the  windward 
streamline  to  a  circumferential  position  of  1  degree 
at  XMAX  (which  will  be  input  as  PHIPD) 

PHIPD  circumferential  position  at  interface  (degrees) 

(for  KBMS  =  0) 

XMAX  axial  position  at  which  the  streamline  shift  is  to 
be  made  (for  KBMS  =  1) (dimensionless ,  ft  or  m) 

ISO  Type  of  Integration  to  be  Performed 

0  for  integration  of  both  boundary  layer  and  streamlines 
1  for  integration  of  streamlines  only. 

Subroutine  PRESS  reads  in  the  number  of  axial  and  circumferential 
pressure  stations  as  well  as  the  pressure  data  (which  must  be  in  the 
form  of  a  pressure  coefficient).  The  pressure  data  is  read  in  one  com¬ 
plete  axial  station  at  a  time.  The  parameters  relevant  to  this  sub¬ 
routine  are  also  shared  with  the  pressure  fitting  routine  and  are  as 
follows : 

NCS  number  of  circumferential  stations 

NAS  number  of  axial  stations 

PHI(j)  array  of  circumferential  pressure  stations 

(NCS  values  to  be  input) (degrees) 

X(j)  array  of  axial  pressure  stations  (NAS  values  to  be 
input) (dimensionless ,  ft  or  m) 
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CP(i,j)  pressure  coefficient  data 

i  =  (1,  NAS),  j  =  (1,  NCS)  . 

The  axial  pressure  derivatives  (at  the  interface  and  body  end)  for  each 
circumferential  plane  are  read  in  two  at  a  time. 

CPX(l,j)  axial  derivatives  at  interface  and  body  end 
CPX(2 , j )  .  =  NC$) 
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2.  Subroutine  BGEQM 

Subroutine  BGEOM  computes  the  geometric  properties  relative  to  the 
body  axes  used  in  the  streamline  and  boundary  layer  calculations.  For 
an  input  axial  position,  subroutine  BGEOM  computes  the  body  radius  and 
its  derivative,  and  the  body  angle  r  and  its  derivative  (see  Figure  2.2). 
A  call  to  subroutine  BGEOM  has  the  form 

CALL  BGEOM  (X,R,DRDX,GM,DGX) 
where  the  input  argument  is 

X  axial  location,  x 

and  the  output  arguments  are 
R  body  radius,  r 

DRDX  dr/dx 

GM  angle  r  =  tan"1  (dr/dx) 

DGX  dr/dx  . 
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3.  Subroutine  BLOTNR 


Subroutine  BLOTNR  evaluates  the  boundary  layer  parameters  for  use 
in  Blottner's  method.  These  parameters  consist  of  C,  the  transformed 
coordinate  along  the  streamline  and  3  the  pressure  gradient  term. 

On  the  spherical  nose,  the  expression  for  f,  is  integrated  by  the 
Runge-Kutta  method  to  yield  f,.+j  and  f„.+1  (see  BOUNDARY  LAYER  METHODS). 
On  the  afterbody,  is  evaluated  by  the  integration  routine  used  in 
the  main  program  and  then  becomes  an  input  to  the  subroutine. 

A  call  to  subroutine  BLOTNR  has  the  form 

CALL  BLOTNR  (PSI ,S,Y,F,DEXI ,EXIH,BETA,DST) 

where  the  input  arguments  are 

PSI  position  on  spherical  cap  in  radians  of  arc  length,  if> 

S  distance  along  streamline  (used  only  on  afterbody),  S 

Y(j)  array  of  dependent  variables  (j  =  1,6) 

F ( j )  array  of  first  derivatives  of  dependent  variables 
{j  -  1,6).  F  -  $ 

DEXI  transformed  step  size  along  streamline,  Af,  =  f.^  +  j  - 
DST  step  size  along  streamline.  As 

and  the  output  arguments  are 

EXIH  transformed  coordinate  at  mid-point  of  interval,  r.  , 

1+2 

BETA  pressure  gradient  parameter,  3 
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4.  Subroutine  FC N 

Subroutine  FCN  computes  the  first  derivative  for  each  of  the  de¬ 
pendent  variables  to  be  used  in  one  of  the  streamline  inteqration 
routines. 

A  call  to  subroutine  FCN  has  the  form 

CALL  FCN  (N,S,Y,F) 
where  the  input  arguments  are 

N  number  of  differential  equations 

S  distance  along  streamline  (independent  variable),  S 

Y(j)  array  of  dependent  variables  (.i  =  1,N) 
and  the  output  argument  is 

F(j)  array  of  first  derivatives  (j  =  1,N),  F  =  . 
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5.  Subroutine  COEFF 


Subroutine  COEFF  calculates  the  coefficients  to  the  appropriate 
finite-difference  boundary -layer  equations.  These  coefficients  are  func¬ 
tions  of  both  the  normal  and  tangential  grid  spacing,  the  pressure 
gradient,  and  the  velocity  profile  at  the  previous  computational  station 
(see  BOUNDARY  LAYER  METHODS). 

A  call  to  subroutine  COEFF  has  the  form 

CALL  cor  I r  ( JMAX ,AM,BM,CM,DM,AS ,CS ,DS ,W,WL,KB) 
where  the  input  arguments  are 

JMAX  number  of  grid  locations  in  normal  direction 

W(j)  matrix  of  present  iterative  values  of  transformed 

velocity  components  (j  =  l,JMAX)(dimensionless) , 

Vl,j 

WL(j)  matrix  of  transformed  velocity  c  nponents  at  last 

integration  station  (j  =  1 ,JMAX)(dimensionless) , 

Ui 

KB  indicator  variable 

1  for  calculation  of  coefficients  of  Blottner's 
similarity  equations  at  stagnation  point 

2  for  calculation  of  coefficients  of  Blottner's  non¬ 
similar  equations 

3  for  calculation  of  coefficients  of  Hall's  equations 
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and  the  output  arguments  are 


AM( j )  ' 

BM(j) 

CM(  j ) 

DM(j)  •  coefficients  of  respective  boundary  layer  equations 
AS(  j ) 

CS 

DS(j)  , 

(all  arrays  are  of  dimension  JMAX) 
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6_. _ Subroutine  INVERT 

Subroutine  INVERT  solves  a  block  tridiagonal  system  of  linear  equa¬ 
tions  using  the  modified  Davis  algorithm.  The  coupled  continuity  and 
momentum  equations  form  such  a  system. 

A  call  to  subroutine  INVERT  has  the  form 

CALL  INVERT  (JMAX,A,B,C,D,AS,CS,DS,W,KB) 
where  the  input  arguments  are 

JMAX  number  of  normal  grid  points 


coefficients  to  respective  boundary  layer  equations 
(all  arrays  are  of  dimension  JMAX) 


KB  indicator  variable 

1  for  Blottner's  method  at  stagnation  point 

2  for  Blottner's  method  at  all  other  points 

3  for  Hal  1 ' s  method 


A  ( j ) 
B(j) 
C(j) 
n(j) 

AS  ( j ) 

cs 

DS(j) 
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7.  Function  KRUNGE 

Function  KRUNGE  is  a  subprogram  which  uses  the  fourth-order  Runge- 
Kutta  method  to  integrate  a  system  of  NDE  first-order,  ordinary  differ¬ 
ential  equations  with  a  variable  step  size.  As  a  criterion  for  varying 
the  computing  interval,  the  differential  equations  are  integrated  over 
an  interval  of  step  size  DSS  first  and  then  over  the  same  interval  with 
two  step  sizes  of  DSS/2.  The  two  solutions  are  then  compared  to  give 
an  estimate  of  the  error  for  each  variable.  If  any  error  is  larger  than 
EPS  =  E  -  04,  these  answers  are  discarded  and  the  comouting  interval  H 
is  halved.  If  all  of  the  error  estimates  are  less  than  EPS,  the  answers 
are  allowed  and  the  integration  process  continues.  In  addition,  the 
step  size  is  either  doubled  or  set  equal  to  DST,  whichever  is  the  smaller, 
for  the  next  integration  cycle. 

The  function  RUNGE  is  used  in  the  main  program  and  has  the  form 
K  =  KRUNGE  (Y ,F,S ,DSS , NDE , DST ,MR) 
where  the  input  arguments  are 

Y(j)  array  of  dependent  variables  to  be  integrated, 

(j  =  1,NDE) 

F(j)  array  of  first  derivatives  of  the  dependent 

variables,  (j  =  I, NOE),  F  =  dY/dS 

S  independent  variable  (distance  along  a  streamline) 
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DSS  integration  step  size 

NDE  number  of  differential  equations 

DST  maximum  integration  steo  size,  AS 

MR  indicator  variable,  MR  =  1  for  the  previous  inte¬ 

gration  interval  to  be  recomputed  with  a  new  step 
size  DSS  determined  in  the  main  program 

and  the  output  arguments  are 

Y(j)  array  of  updated  dependent  variables  (i  =  1,NDC) 

F(j)  array  of  updated  derivatives  (,i  =  1 , NDE ) 

K  indicator  variable 

0  implies  completion  of  integration  cycle 

1  implies  the  integration  cycle  has  not  been 
completed 

2  implies  the  step  size  has  been  reduced  to  a 
value  below  E-08  . 
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8.  Subroutine  MI LNES 

Subroutine  MILNES  uses  the  fourth-order  predictor-corrector  method 
of  Milnes  to  numerically  integrate  a  system  of  NDE  first-order,  ordinary 
differential  equations.  Since  this  method  is  not  self-starting,  it 
must  be  used  in  conjunction  with  an  alternate  method  (such  as  the  Runge- 
Kutta  or  Gear  method)  to  generate  the  starting  values.  The  advantage 
of  a  method  of  this  type  is  that  the  computational  work  is  keDt  to  a 
minimum  between  integration  steps.  This,  however,  is  achieved  at  the 
price  of  accuracy  since  the  step  size  is  held  constant  over  the  entire 
interval  regardless  of  the  error  introduced. 

A  call  to  subroutine  MILNES  has  the  form 

CALL  MILNES  (Y,F,PCM,NDE,DST,S) 
where  the  input  arguments  are 

Y(j)  array  of  dependent  variables  (j  =  1,NDE) 

E(j)  array  of  first  derivatives  of  dependent  variables 
(j  =  1  ,N0E) ,  F  =  jj-*- 

PCM  temporary  storage  of  both  the  dependent  variables  and 

their  first  derivatives  at  four  previous  stations 

NOE  number  of  differential  equations 

l)ST  step  size  along  streamline,  AS 
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S  streamline  distance  (independent  variable), 

S 

and  the  output  arguments  are 

Y(j)  updated  array  of  dependent  variables  at  next 
integration  step,  (j  =  1,NDE) 

F(j)  updated  array  of  first  derivatives,  (j  =  1 ,NDE) 

S  updated  value  of  independent  variable,  S  . 
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9 . _ Subroutine  DGEAR 

Subroutine  DGEAR  integrates  a  system  of  first-order  differential 
equations  using  the  backward  differentiation  formulas  of  Gear  (Ref.  12). 
This  technique  is  particularly  well  suited  to  situations  in  which  the 
system  of  differential  equations  may  be  classified  as  stiff.  In  these 
types  of  applications,  other  techniques  would  be  apt  to  decrease  the 
integration  step  size  to  prohibitively  small  values  in  an  attempt  to 
satisfy  the  allowable  error  tolerance.  Gear's  method,  however,  has  the 
property  of  "stiff  stability"  which  effectively  removes  the  limitations 
on  the  step  size.  The  integration  step  size  is  adjusted  in  the  routine 
so  as  to  satisfy  the  error  tolerance  specified  by  the  user.  The  tech¬ 
nique  used  is  similar  to  that  employed  in  the  Runge-Kutta  method  de¬ 
scribed  in  Function  KRUNGE.  The  method  also  necessitates  that  (in 
general)  a  nonlinear  system  of  algebraic  equations  be  solved  at  each 
step  of  the  integration.  To  solve  these  equations,  the  integration 
package  has  the  option  of  employing  a  variety  of  iterative  schemes. 

A  call  to  subroutine  DGEAR  has  the  form 

CALL  DGEAR  (NDE ,FCN,FCNJ ,S,HG,Y , SEND, TOL ,METH , MITER , INDEX ,IWK,WK,IER) 
where  the  input  arguments  are 

NDE  number  of  differential  equations  to  be  integrated 

FCN  subroutine  to  evaluate  first  derivatives  of 

differential  equations 


/ 
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FCNJ 


S 

HG 

Y(j) 

SEND 

TOL 

METH 

MITER 


subroutine  to  evaluate  the  Jacobian  of  the  system 
of  differential  equations  --  this  parameter  may  or 
may  not  be  specified  depending  on  other  quantities 
specified  in  the  argument 

distance  along  streamline,  S 

integration  step  size  along  streamline 

array  of  dependent  variables  at  present  station 
(j  =  UNDE) 

value  of  independent  variable  at  which  dependent 
variables  are  desired  ,  S  +  AS 

maximum  error  tolerance  allowed  between  integra¬ 
tion  steps 

Basic  Integration  Method 

1  for  use  of  Adam's  method 

2  for  use  of  Gear's  method 

Iteration  Method 

0  for  functional  iteration,  internal  calculation 
of  the  Jacobian 

1  for  chord  method,  Jacobian  is  supplied  externally 

2  for  chord  method,  internal  calculation  of  Jacobian 

3  for  chord  method,  diagonal  approximation  of 
Jacobian  is  made  internally 
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1 

! 


INDEX  1  for  first  call  to  subroutine 

0  for  remaining  calls 

IWK  work  vector  of  length  NDE 

WK  work  vector  of  length  13  x  NDE 

IER  error  parameter 

33  for  error  test  not  satisfied  due  to  too  low 
an  error  tolerance 

66  for  error  test  was  satisfied  only  after  HG 
was  reduced 

132  error  test  failed  after  HG  was  decreased  to 
lower  limit  . 
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10,  $ ubroutine  PRESS 

Subroutine  PRESS  reads  in  the  pressure  data  used  in  the  pressure 
fitting  technique.  The  input  pressures  must  be  in  the  form  of  a  pres¬ 
sure  coefficient.  The  data  read  into  the  subroutine  consist  of  the 
following: 

NCS  number  of  circumferential  pressure  stations 

NAS  number  of  axial  pressure  stations 

PHI ( j )  array  of  circumferential  stations  (degrees) 

(j  =  1,NCS) 

X ( i )  array  of  axial  stations  (i  =  1,NAS) 

CP(i,j)  pressure  coefficient  data  (i  =  1,NAS,  j  =  1,NCS) 

CPX ( i , j )  axial  pressure  derivatives  at  body  interface  and 
end,  (i  =  1,2,  j  =  1,NCS)  . 
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11.  Subroutine  SPHCAP 

Subroutine  SPHCAP  computes  the  inviscid  flow  prooerties  on  the 
spherical  nose  of  the  body  to  be  used  in  the  boundary  layer  integra¬ 
tion.  On  the  initial  call  to  the  subroutine,  a  check  in  continuity  be¬ 
tween  the  pressure  and  its  axial  derivative  across  the  interface  is  first 
made.  If  the  values  of  the  parameters  should  vary  by  more  than  0.003  or  4.0,  re¬ 
spectively,  the  quadratic  spline  technique  is  extended  to  the  stagnation 
point.  Otherwise,  the  series  expression  for  the  pressure  coefficient 
will  be  used.  This  in  turn  will  yield  the  dimensionless  fluid  velocity 
which  is  used  in  the  subsequent  boundary  layer  calculations.  Additionally, 
the  derivative  of  the  pressure  coefficient  with  respect  to  ip  is  calculated 
for  use  in  Blottner's  boundary  layer  method. 

A  call  to  subroutine  SPHCAP  has  the  form 

CALL  SPHCAP  (P ,CP ,UE ,DCPSI ) 
where  the  input  argument  is 

P  position  on  nose  in  radians  of  arc  length  form 

stagnation  point,  ip 

and  the  output  arguments  are 

CP  pressure  coefficient,  Cp 

IJF  fluid  velocity  at  edge  of  boundary  layer 

(dimensionless),  ue 
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OCRS  I  derivative  of  pressure  coefficient  with  respect 
3C 

*•  i  ■ 

an  additional  input  to  the  subroutine  is 

COE(i)  coefficients  used  in  the  series  expression  for  CP 
(i  =  1,10)  and  is  specified  in  a  data  statement 
within  the  subroutine  (see  SURFACE  PRESSURE  DISTRI¬ 
BUTION). 
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12.  Subroutine  MIDPTS 


Subroutine  MIDPTS  calculates  the  dependent  variable  at  the  mid- 

♦  point  at  each  of  the  subintervals  defined  within  the  parabolic  spline 
technique  (see  SURFACE  PRESSURE  DISTRIBUTION).  This  amounts  to  solving 
a  tridiagonal  system  of  linear  algebraic  equations  by  the  LU  decomposi¬ 
tion  method. 

A  call  to  subroutine  MIDPTS  has  the  form 

CALL  MIDPTS  (AA,BB,CC,D,YM,N) 
where  the  input  arguments  are 

arrays  of  coefficients  of  tridiagonal  system 
(i  -  1,  N  -  1) 

number  of  discrete  data  points  along  interval  to  be 
spline  fit 

>  and  the  output  argument  is 

•  YM( i )  array  of  dependent  variables  at  midpoints 

(i  =  1,  N  -1),  Y  . 


AA(  i ) 
BB(i ) 
CC(i) 
D(i)  . 

N 
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13.  Subroutine  INVISD 


Subroutine  INVISD  computes  the  properties  of  the  inviscid  flow 
field  at  any  point  on  the  body  surface  after  the  interface.  These  pro¬ 
perties  consist  of  the  pressure  coefficient*  the  local  fluid  velocity 
at  the  edge  of  the  boundary  layer,  the  first  and  second  ci rcumferential 
pressure  coefficient  derivatives,  the  axial  pressure  coeffic.ent  deriva¬ 
tive  and  the  mixed  derivative.  These  parameters  result  from  a  quadratic 
spline  fit  to  the  input  pressure  data  (see  SURFACT  PRFSSURE  DISTRIBUTION). 
The  method  employed  is  limited  in  application  to  situations  in  which 

the  number  and  position  of  the  circumferential  stations  do  not  vary 

between  axial  stations. 

A  call  to  subroutine  INVISD  has  the  form 

CALL  INVISD  (XX,PPH,UE,DCPX,DCPPH,CPC,DCPXP,02PPH) 
where  the  input  arguments  are 
XX  axial  position,  x 

PPH  circumferential  position  (radians),  <|< 

and  the  output  arguments  are 

UE  fluid  velocity  at  edge  of  boundary  layer  (dimen¬ 

sionless),  ug 


DCPX 


3C 

axial  derivative  of  pressure  coefficient, 


DCPPH 


CPC 

DCPXP 

D2PPH 

other  inputs 
NAS 

NCS 
X  ( i ) 
PHI(j) 

CP(U) 

CPX(i,j 


1 

circumferential  derivative  of  pressure  coefficient, 


pressure  coefficient,  Cp 

mixed  axial -circumferential  second  derivative  of 

3ZC 

pressure  coefficient, 

second  ci rcumferential  derivative  of  pressure 

3?C 

coefficient, 

to  the  subroutine  consist  of 


number  of  axial  pressure  stations 

number  of  circumferential  pressure  stations 

array  of  axial  pressure  stations  (i  =  1,NAS) 

array  of  circumferential  pressure  stations 
(degrees),  (j  =  1,  NCS) 


pressure  coefficient  data  (i  =  1,NAS  ,  j  =  1,NCS) 

)  axial  pressure  derivatives  at  body  interface  and 
end  (i  =  1,2  ,  j  =  1,NCS)  . 
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14.  Subroutine  STAG N 


Subroutine  STAGN  locates  the  stagnation  point  on  a  given  configura¬ 
tion.  For  the  case  of  the  sphere-ogive-cylinder  geometry,  this  amounts 
to  calculating  the  Newtonian  stagnation  given  an  effective  angle  of 
attack.  For  configurations  having  analytical  pressure  distributions, 
the  stagnation  point  may  be  calculated  from  analytical  expressions  which 
must  be  supplied  by  the  user  (as  was  the  case  of  the  ellipsoid  of 
revolution) . 

A  call  to  subroutine  STAGN  has  the  form 

CALL  STAGN  (ALP.XO.XNOSE ,RPFR) 
where  the  input  arguments  are 

ALP  effective  angle  of  attack  (radians), 

XNOSE  distance  from  oriqin  of  body  axes  to  body  nose 
RPER  radius  of  spherical  cap,  Rper 
and  the  output  argument  is 

XO  axial  location  of  stagnation  point 
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15.  Listing  ot  Main  Pcogras 


♦ 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

e 

c 

c 

c 

c 

c 

c 

c 

c 

c 

C 


DIMENSION  AM  (50)  ,  Oil  (50)  ,CM  (50)  ,  Ult  (50) 

DI8E  KS1CN  US  (50 )  ,  U(2,50),  UL(2,50),  U0(40) 
DIMENSION  AS  (50)  ,DS(50) 

DIMENSION  Y (6) ,F (6)  ,  PCM (6, 4,  i)  ,IRK  (6)  ,  «K  (UO) 

EX  IE  UK  A L  FCR,  FCN  J 
CC6.1CN  /HALL/DSTN,DY,UE 
CO anC  K  /BLCl/CKTA,  OEXI ,  LXIll,  BIPB 
CC lines  /SC  ALL/1IL,  ilLii ,  H 

CCF.MCK  /INTFc/XINT.PHIP.S  AL,CAL,PSIMAX,  PSIU,RP£E 
CCBJ1CN  /CPSTA/NAS.NCS 

COMMON  /OUTPT/&,CF,DCPX,DOPPH,  DOP XP , D2P PH , GH , DGX 


DESCBIFTICN  OF  INPUT  PAHAMETEBS: 


BPEB  E  A 0 1  US  OF  SPUERICAL  CAP 

XNOSE  EIS1ANCE  FRCfl  ORIGIN  CF  BODY  AXES  Tu  NOSE  OF  BODY 
BL  TOTAL  LOO Y  LENGTH 

XI  NT  AXIAL  ICC  A  TICK  CF  S  PHEN  E-.A  FT  Fn  BCD  Y  INTEBFACE 

DST  INTEGRATION  STEP  SIZE  ALONG  SIBEABLI ME 

DPSI  STEP  SIZE  ON  SPHERICAL  CAP  ( CtGREES  OF  ABC  LENGTH) 

ALP C  EFFECTIVE  ANGLE  OF  ATTACK  (DEGREES) 

NBS  C  FOG  HALL'S  BOUNDARY  LAYEfi  EET HOD 

1  FCB  ELOTXNEB'S  METHOD 

NT  EDGE  TEST:  0  FOB  POINT  ADDED  AT  EDGE  HHEN  VELOCITY 

AT  POINT 

NEXT  TO  EDGE  IS  LESS  THAN  99.95  PEBCENT  OF 
POINT  AT  EDGE. 

1  FOB  POINT  ADDED  « HEN  EDGE  SKIN  FRICTION 
COEFFICIENT  IS  ABOVE  0.0005. 

NO  CONVERGENCE  TEST:  0  FOE  CONVERGENCE  BASED  ON  CFBEX 

CHANGING  BY  LESS  THAN  0.5  PEBCENT 
BETUEEN  SUCCESSIVE  I1EBAUONS. 

1  EACH  PCI NT  IN  PROFILE  IS  CHANGING  BY 
LESS  THAN  0.1  PEBCENT. 

nor  METHOD  USED  10  INTEGBA1E  STREAMLINES: 

1  FCB  HILNE'S  PREDICTOR  CQKRECTG8  METHOD 

2  tCii  tiUNGE-KUTTA  METHOD 

3  FCR  OtAh'S  METHOD  (STIFF  SYSTEM) 

MAXS  NUMB  Eu  CF  STATIONS  TO  b£  COMPUTED 

K P=  fi  FCB  FuCf ILLS  PRINTED  EVERY  NTH  STATION 

K PH  C  FCR  NO  INTERATION  PBOFILES  r BINTED 

1  FOB  ITERATION  PROFILES  PRINTED 

KPO  0  FCR  ORDINARY  PRINT  OUT 

1  FCR  ALDxl ICN AL  PRINT  OUT 

KBM  NUMBER  CF  STREAMLINLS  TO  TR  ACL 

KDMS  C  FCR  ORDINARY  BUN 

1  FCR  STREAMLINE  oUItT  FROM  EUI-0  TO  PHI-1.0  AT  XMAX 
PHI  PD  CIRCUMFERENTIAL  ANGLE  AT  INTERFACE  (KbMS-U)  (DEGREES) 
XMAX  PCSITaO i  CN  WINDWARD  STREAMLINE  U  HER  E  STREAMLINE 

SHIFT  FROM  WINDWARD  PLANE  IS  MAOL  (KEMS*1) 
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ono  noonnn  n n 


C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


RE  AC  (1,2)  FPE B,  mo St,  3L,XINT,DST,  CBS  ID,  ALPD 
2  FORMAT  (4  (IX, F  10.5) 3 (IX, P 10.  5) ) 

h£  AC  (1,4)  NBS#NT,NC,HOT,MAXS,KP,KPH,KPO,KBn,KBeS,ISO 
4  FORM  AT  (11  (  1  X , 13) ) 

JMAX=4C 
£T  AE-4 • 5 
YECGc*C.  15 
CALL  CRESS 
R0E  =  5-»KES 

if  USC.  EQ.  1)  GO  TO  36 
46  *(1,1)=0.0 

W  (1,JMAX)  *  1.C 
W  (2,  1)=0.0 
•  L<1, 1)30.0 
WL  ( 1 , J P  AX)  *1.0 
*LU,1)-=0.C 
DE1A=£T AE/  (JflAX-  1) 

I  i  E  h  *  0 
CEbEXL’I. C 
C Y*  YE  CO  E/  (JHAX-  1) 

D YS»CY 
DSTS*CST 


ISO  1  FCE  INTEGRATION  Of  STREAMLINES  ChILY 

C  FCE  ADDITIONAL  INTEGRATION  OF  bCUNDAFY  LAYER 


STREAMLINE  INTEGRATION  : 

S  «  STREAMLINE  CISTANCE  (INDEPENDENT  VARIABLE) 
1(1)  AXIAL  ECSII10N 
1(2)  CIECUHf EREN'TIAL  POSITION 
Y (3)  STREAMLINE  ANGLE 

Y  (4)  L  (T££TA)/D  (PHI) 

Y  (5)  LOG  (D  (Eiix)/0  (BETA)  ) 

Y (b)  EXL  (FOR  BLOTTNER * S  METHOD  ONLY) 


P1=ACCS  (-1.) 
DGR=PI/180. 


GUESS  INITIAL  PROFILE 
DO  10  0*2,  UMAX 
a  (i,J)*i.C 
h  L  ( 1 , J ) >1.0 
■  S (0) 3  1. o 

h  U,J)  *N  (2,J-1)  -0.5*DETA*  (U  (1,J)+N  (1,  J-1)  ) 
10  CENTIME 


Kb*  1 


112 


30  CONI IN  IE 

CALL  CC2FF  (oflAX, Aa,brt,ca, DB,  AS  ,CS  ,  CS  ,  V. ,  WL,  K3) 

CALL  INVERT  (JMAX,A(i,bfl,CH,Dfl,AS,C'S,i>S,li,Kb) 

1090  F08BS1 (/, 10 (2X,F8.  6) ) 

ITEF  =  I1EM  1 
KB  =  2 

IF  (AC. EL.O)  CO  TO  32 
DO  34  J  =  2,  J h AX 

IF  (M'S  (  (US  (J)  -y  (1  ,J)  )/U  11  ,J) )  .GE.  0.001)  KD*1 
WL(1,J) =  W  (1,0) 

34  US (J) = h  (1  ,  J) 

IF  (K £.  K.  1)  GO  10  30 
32  FP  =  (2.* Vi  (1,2) -0.5*4  (1,3)  )/DETA 

CEhEX  =  i.*£cH  U.) *FP 
CFREXC*CF  RE* 

IF  (NO. Ft.  1)  00  TO  36 

IF (  (AES  (CFBti-CFNEXL) /CFfeEXL)  . GT.0.U05)  KB=1 

CF  RE  XLSCF  BEX 

IF  (KB. EC.  1)  00  10  30 

FPE*“ (2.  *k (1,JMAX-1)-0.5*» (1.JMAX-2) -1. 5*H (I.JKAX) ) /DETA 
CFBEXE*2« *S*ET (2. ) *FPE 
IF (CFBEXE.LE.O. 002)  GO  TO  36 
ETAE=E1AE4CEIA 
GO  XC  4b 
36  rfBiTE  (2 ,5| 

5  FOBHA'1  (//,20a,'  INPUT  PARAMETERS:  •  ,/) 

.RITE  (2,3)  EPER,XNOSE,BL,XI»T,DST,£i'Si.D,ALPD,aAXS,NT,NC,BOT, 
*N3S,KP,KPH,kE0, NCS.NAS 

3  FORflAT  (/,23l,  'F.PEB  =  '  ,bX,F0.  5, 5X,  'X  (NOSE)  =  •  .  5X,  F8.  5, /,  23X, 

•  * 3C  D  Y  lENGld  =  ', IX, F8. 5, 5X, 'X  (INTERFACE)  =  • ,F8. 5,/,2JX, 

**061  =  • ,E*,ffa.S,5X, 'DELIA  PS1  =  ' , 3X , F8. 5, /, 2JX, 

•  'MiG,  ATTACK  *  '  ,  1  X,  F8. 5,  5X,  •  NC  .  STATIONS  *  ',I8,/,23X, 

* ' E  CO  E  1  ESI  =  ' , 3X,I2, 3X, 'CCNV.  1EXT  » * , 2X,I2, 3a, 1 INIEG.  BET 8  =  • 
*,lX,Ii,/,23X,‘D.L.  METHOD  =  • , U, 12, 3X, » NC.  PBOFILES  “  *,I2,3X, 

•  •HER.  PBCfiLES  **  ,I3,/,23X,*jia1RA  P.O.  *  *,21,12,3  X, 

*/  ,23  X, ' NO.  AXIAL  PRES5UBE  STATIONS  *  ',I3,/,23X, 

••NC.  CIRCUrtt EbENTIAL  PRESSURE  STATIONS  *  *,13,//) 

IF (NCS. Ew.O. 04.  KBB. EC* 1* Oh.ISO. Eg.  1)  GO  TO  45 
OC  4?  J- 1 ,JH AX 
»0  (J) =,  (1  ,J) 

45  CCaIIKUE 
OElAb-LETA 
ALL  =  AIFC*CGl« 

C  A I  =  C  C  S  (Alf) 

SAl*SIN  (ALE) 

yc=o.c 

G  fl  C  =  9  C  .  *CGR 

CALL  S1A0N  (AIP,XO,XNOSE,hPEK) 

IF  (XC.NE.O.U)  CALL  UG  EOfl  (XO  ,  YO  ,  DR  C  X  ,  GBO  ,  DG  X) 

CGC  =  CC £  (GBO) 
jGC *£I  N  (GHCJ 
wRITE  (2,7)  aC, Yu 

7  FOpHAI  (/ , 2 JX , • SI Aw.  POINT,  *0  = ' , Ptt, 5 , / , 23 X , • STAG.  POINT,  10  =  ', 
•Fd . 5 ,/ ) 

IKKF.Fw.1l  *f.ITE  (3,  1090)  («  ( 1  ,J)  ,  J«  1 ,  JBAX) 

■  RITE  (3,7)  C i REX , ITER 

9  FGhilAT  (/,  1CA, 'STAGNATION  POINT  CFhFX  *  ',F9.6,2X, 

* ' CC  N  V  £  FUE  NC  t  IN  ' , 13, IX, ' ITERATIONS' ,/) 

P3  R  =  ACCS'  (  (6CE«-XINT-fXNoSE)  /h PEE)  -  ALP 
IF  (FSIC.Gl.LGR)  PSIC=CGR 
Fj1CC*ESIC/0Gp 
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n  n 


33 


UCB=SAb  (PSIC) 

D£ES*LCfWbE£fi 

CALL  BGECH  (XINT , hi , DPI ,GNI, DU1 ) 

SGI=SIb  (GRI) 

CG1=CCS  (Gfll) 

■  SITE  (3,33)  D0K,PS10D, DESID 

FORMAT  {/,  10Xf 'DtPS/KPEH  «  • , F14. 6,/, 1 Ox, • PSI  (EPS) 
♦•DESI  =  '  ,  f  5. 2,/) 


=• ,F7.3 


KB£=  1 

8  REAL  (1,6)  EHIPO 
o  FCbMAT  (F10.5) 

If  (KERS.Evi.O)  GO  TO  12 
XM AX=EE1PD 
PUIEU  =  C .  0 

12  CCb'IIKOE 
JMAX*40 

IF  (KES.EQ.  1.UE.3BS.EQ.0)  GO  TO  13 
DC  13  J=1,JHAX 
WL  (1 ,J) =WC  (J) 

13  CCbTlbLE 
LK-0 
HL»0. 

Ki=0 
NHE=  0 

c*i»cmxc 

S=6EEE*ESIC 
SI* G  •  C 

PHlE=EHIPO*UGh 

DENS  vs  ( 1  •  +  T  AN  (GHI )  **2)  *SIH  (PHIP)  •  *2*TAN  (ALP)  **2 
♦-HI.  +1AMGBI)  *CQS  (PHIP)  *TAN  (ALP))  **2 
BETA3  ACCS  (  (CCS  (PHIP)  +T  AN  (GHI) *TAH  (ALP)  )  /S^iiT  (DENSy)  ) 
3EIA£=EETA/DGfi 

PSIP.  I  X=  ACCS  (CAL*SGl+SAL*CGI*COS  (PHIP)  ) 

PS I A Cs ESIN AX/ OGB 

DES1=CESIC*CGH 

FSI=P£IC 

IF  (ESl.GE.ESIMAX)  PSI-PSI3AX 
NiiPS*  (PS1MAX-PS10)  /DPS I 
unlit  (3,31)  EHIPD,BETAE,PSIMD 

51  ECRHAl  (/, 1CX,'PHI  (INTERFACE)  «  • , F6. 2 ,/ , 10 X, • BET  A  =  ',F7. 
*10X,*ESI<R4X)  «  * , F 8 . 4 , // ) 

CALL  SEBCAE  (0. 0  ,CP  ,  U1S,  DCPSI) 

If  (ISC.  tv).  0)  GO  TO  53 
H=bEER*SIN  (ES1MAX) 
riL*H 

S*H  EB*£SIHAX 
GO  TO  42 
53  CCNTINOE 

4 R IT  E  (3,35) 

35  FORMAT  (//,  10 X, • SPHERIC AL  CAP  RESULTS:  *  ,/) 

■  RITE  (3,52) 

52  fCRH&l (tX,'ESI* , 13X,'CP* , 1 4X, » UE/V lb • , 10X, « CFR lx • , 1 IX , • H‘ 
*•3',  15X,' I1ER', 12X,' CFREX  (EDGE)  • ,/) 

DP5=CESI*fEt6 

IF  |.»eS.Ey.  1)  GO  TO  39 

CALL  SEHCAE  (ESI,CP,UE, CCPSI) 

DC  oO  J-2, JMAX 
WL  ( 1 , J ) *0 . C 
w  (2 ,  J)  *  0.  0 


,/»10X, 


3,/. 


,  15X, 
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<i  (1,J)  * UE 

60  CCMISUE 
GO  TC  4  1 

39  CALL  SEHCAP  (PSI,CP,0E,0CPSI) 

41  H=EEER*SIN  (PS1) 

Nb£=NHl4l 
GC  10  6C0 

40  PSIC  =  P  £  1/ DGR 

WBl'IE  (3,37)  PS1D,CP,UE,CE'BEX,H,S,  1TLR, Cf B£XE 

3  7  F0BBA1  (/,  2X,  C  (3X,  F  12.  5)  ,  1  IX, 12, 3X , FI  2 . 5) 

IP  (Kl. .Eg.  1)  GO  TO  42 

IF  (SHE. II. NUPS)  GO  TO  38 

3=RFER*ESISAX 

PS1=ES1MAX 

KL=  1 

GO  TC  39 
38  S=S-)  CPF 

ESI=P£I4DESI 
GO  TC  39 
C 
C 
C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

42  LK=1 

HG=0. COOO  1 
TCI*  C. COOO  1 
Milt  b*  3 
INC£X=  1 
M  E  T  H  *  2 

»S1'1E  (3,1  C 40) 

1040  £  0  DM  A1  (//  ,  6  X  ,  'X* , 1 5X, *  PHI ' , 13X,' THETA', 11X,'S  ',14X, 
*'H»,1SX,'UE /VIS' , 10X,'DPDS,,!2Xf,CFDEL’) 

IF  (KIC. EQ.  1)  4RITE  (3,  1042) 

10  42  FoiulAl  (6X  ,  '  CE  '  ,  1 4X, ' DCP/U  X'  10X  ,  '  DCF/DPHI  •  ,  8X  ,  '  DCPXE*  ,111, 

A'Dif  Pli  •  ,  1 1X,  •£• , 15X, ' ITEB  *,  10X,  ’CFDEL  (EDGE) ' 

AbX,  •  C  (TfcElA) /LPHI',3X,'LGu(DPHl/DEETA) •) 

1050  FC8MAI |//,2a,4(3X,F12.5) ,  )X,E14.6,3(3X,F12.5)) 

10  52  F08MA1  (2X  , 6  (JX,  F U.  5)  ,  3X  ,  1 10, 5X ,  P 12.  5,/,2X,  2  (3X,  FI 2.  5)  ) 

C 

C 

C 

C  UEG1N  1C  IbACE  SIBcAMLINE  OVER  DOCK,  AFTER  EACH  INTEGRATION  STEP 

C  COM  PL 1 E  IDE  ECUHLARY  LAYER  USING  EITHER  ULCTTNE6S  OR  HALLS  METHOD 

C 

C 

C 

C 

C  INIT1A1  VALUES  AT  INTERFACE 
LS1=C£1S 
DY=CY£ 
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ooodo  non 


L'ET  A  =  C  El  AS 
C 

Y  (1) =  X1 NX 
3PSI*S1N  (ESittAX) 

CPSI=CCS(ESIBAX) 

Y  (2)  *EBIP 

Y  (3)  =ACCS  (  (CAi*CGI-SAL*SGI*CO:i  (Y  (2)  j  )  /SFSI) 

If  (Y  (2)  .If. 0.0)  Y  (3)  =0.0 

CPSIEE  =  (SAI*CG1*SIM {Y  (2) ) /SPSi) 

If  (Y  (2)  .  L£.  0.0)  GO  TO  43 

Y  (4)  =  (SAL*SGI*SI»  (I (2)  J-CQS  (Y (3) > *CPS1* DPSIDP) / 

A  (-SPSI*SIN(Y(3)|  ) 

43  I F  ( Y  (1)  .LI.  0.0)  Y  (4}  =SQiiI  (-i»AL*SGI/aPSI+C?SI  *SAL*CGI/SPSI**2) 
Y (5)  *  A  ICG (StfcB*SPSI/ (BI*CGS  (Y  (3) ) ) ) 

IF(NES.EC.O)  GO  TO  44 
CALI  ECU (HDE,S,Y,F) 

CAU  LlCTN&(FSX,S,Y,F,EEXl,EXlti,I>lPH,DSI) 

44  K-0 

LP*1 
J  J*3 
LSf *DST 
KL  =  0 

IF (ISC. EC. 1)  GO  TO  130 
IF  (HOT .HE. 1)  GO  TO  550 
DC  3C7  J»1,  ACt 
PCB  (J,4, 1)  *Y  (J) 

307  PCB  (J, 4,2) *0. 


550  3EHOS+CST 

IF (KEHE.EC.O)  GO  TO  560 

IF (Y ( 1) . LT . XflAX)  GO  TO  5oG 

Y  (2)  =£GR 

Y (3)  =Y  (4)  ♦  CGfi 

INCEX=  1 

I1G  =  C .  CGOO  1 

CAU  FCMHCE,S.Y,F) 

XMAX=BL 
560  CCU1ACE 
KHEaNHE+1 

GO  TC  (305,310,315) ,KOT 
305  IF  (JJ.GE.  1)  GO  TO  315 
309  CALL  MILHIS  ( Y ,F , PC3 ,HDE, U3T,S) 
GC  TC  555 

308  DO  580  J«1,N£E 

PCB  (J,JJ,  1)  *Y  (J) 

580  pen  (J,  JJ,2)  af  (J) 

JJ=JJ-1 
GC  TC  555 


310  K  =  Kt>UNGE(Y,F,S,CS3,NDE,bST,Mft) 

IF (K. EC*  3)  GC  TO  142 
CALL  tCA (BC£,S, Y , F) 

IF  lH.Es.-T)  GO  TO  310 

IF (S.GI. SINC. Oi. XL. ty. 1)  GO  TO  311 

SF»S 
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n  n  n  n  n  n  n  non 


00  10  210 

311  If  (Kl.Eg.  1)  GO  TO  312 
NR=  1 

D££=D££*  (SEiiL-SF)  /  (S-SP) 
KL=  1 

00  10  310 
312  KL  =  0 

00  10  555 


315  CftLL  E  G  E  A  h  (NEE,  FCN, FCHJ,5 , HO, I ,  SEkC, T0L,HE7H,MITEB , INDEX , IBK , 

CAII  i CN (N:i(S, i, F) 

IF (rtOl. EQ. 1)  Go  10  308 


555  HaB*CCS  (X  (3)  )  *EXP  (I  (5)  ) 

IF  (ISO .  EQ. 1)  GO  TO  130 

600  LP=  1 

If(KP.KE.O)  LI=MOD  (NHP,XP) 


IF  (N8S.EQ.0)  GO  TC  85 
300  11EE=C 
KB  =  2 

CAII  EICTNo tPSI,S,X,F, DUX  I, EXIH , BIPU , DSX) 

C  liBITE  (3,5CC0)  I  (6)  ,  EIIH ,  DE  XI ,  UlPti 

5000  F0BHA1  (/,4X,  'EXI  =  •  ,  E14.  fa,/,  4  X,  '  EXI H  *  • ,  El  4. 6  ,/,  4X,  •  DEXI  =  •, 
*E14.6,/,4X,'B£TA  *• ,PlU.b,//) 

361  ITEB=0 

360  C AIL  CCEFF(JtfAX,AM(BflfCNfDfl*ASfCSrDS,tf#HLrKB) 

CALL  IKVEE1 ( JflAX, AH, OH ,CH , Dfl , AS , CS , OS , N ,KB) 

IF  (LP.  EL.  C.  ANC.KPH.EQ.  1)  *RI1 E  (3, 1 090)  (H  ( 1 ,  J)  ,  J=  1 ,  JM  AX) 
I'IEB=HEB-f  1 

IF  (ITFE.GE.40)  00  TO  135 
IF  (tiC.  EC.  C)  OC  TC  340 
KU=  3 

DO  25C  J=i,  JbAX 

IF  (ADS  (  (KS  (J) -4  (1  ,J) ) /« (1,J) )  .GE. 0.001)  KB*2 
350  <fS  (J)  =fc  (1  ,J) 

IF  (K2,  tv.  2)  GC  TO  360 
340  FE-(2.*b(1,^)-0.5*b(1,3)) /DEI A 

CFB£X=£gRl  (2.)*tP*SQ&T(3L/!f  (b)  )  *U£**2*H 

IF  (NC.EQ.  1)  GO  TO  370 

KB*J 

IF  (  (ACS  (CF£EXi.~Cf'l>£X)  /CFKEXL)  .OT.  0,005)  KB=2 

CFEEXI=CFFEX 

IF  (Kli.  EC.  1)  OC  TO  360 

370  FP£  =  -(2.*b  (1,j;iAX-1)-0.5*N(1,JHAX-2)-1.5*U  (1, J3AX) ) /DETA 
CffiEXE=CFFEi*F«'E/FP 
i. F  ( til .  EC.  C )  oO  10  375 
IF  (CFFEXE.Gi.  0.0005)  GO  TC  300 
00  10  118 

375  IF  (  («  (  1.JJ1AX-1J  /«  (  1.JMAX)  )  .  LI.  0.  9V95)  00  10  380 
GC  TC  1  18 
380  JHAX=JFAX(1 

IF  (J3AX.G1  .  V) 0 )  00  TO  362 
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tL  |1,0HAX)*(.l{1,GHAX-1) 

DC  39C  0=2, UMAX 
390  a  (1,J)  =  *L  ( 1,J) 

£T  AE  =  £  1  A£4  CrT  A 
UHIIE  (3,1030) 

KB*2 
C  .  L£  = 0 

GO  10  361 
3b2  OH AX  =26 

DO  3c3  J*1,25 
GX=**0-1 

363  WL  (1 ,0) =Hl (1 ,0X) 

«L  (1,26) =UL (1,50) 

El AE=E1 AE4DE1 A 
DETA=£ETA*2. 

DPSI=CESI«2.  0 

cs:=osi*2.o 

EPE=CEE*2. C 

WHI1 E  ( 3 ,  1 19)  OM AX ,  DETA 
DO  364  8=2,JMAX 

364  a (1,0) *iii (1,0) 

KB  =  2 

GC  1C  361 
C 
C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

C  APPLY  HAIL'S  !h  AHS FGPfl ATIOfl 

C  D(X)=C(X/l),  D(Y)=D(SCKT(RE)  *Y/L)  ,  U=U/VIN,  UE=UE/V£N 
85  KB*3 

OSTN  = (6-Sl) /OL 
SN=S/E  L 

3H*SL4 <S-Sl)/2.0 
SHA=SH/BL 
8LH=  (H-fHD/2. 

IF  (*L(1,JHAX) . IE. 0 . )  GO  TO  88 
DC  88  0*2,  Jit  AX 
H ( 1,0)  =  « (1,0)  *UE/«L(1,0HAX) 

88  CC  All  A  CE 
86  1 1  £  h  =  U 

97  CALL  CCfcFE  <0  MAX  ,  AH,  1>M  ,CM  ,  DH,  AS,  CS  ,  CO ,  k  ,  WL,  KB) 

CALL  INVEST  (OMA  X,  AH,  DM  ,CM ,  Drt  ,  AS  ,CS  ,  OS  ,  ii  ,  KO) 

IF  (LE. EC. C. AND. KPu.EL.  1)  WE1TE  ( J,  1090)  (k ( 1 , J) , 0* 1 , OH AX) 
II  Ef  *  11  ER4  1 

IF (1TEP.GE.40)  GO  TO  135 
KB  =  4 

IF  (AC,  Ek.  0)  GO  TO  196 
CO  18*  0*2,  J  HA  X 

IF  (AES  ((WS(O)-K  (1,J))/W(1,0)»  .GE.  0.001)  K0  =  ) 

182  kS|J)  *k  (1,0) 
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IE  (K2. EC- 3)  GO  TO  97 
196  FP=(2.*4(1*X)-0.5*W(1,3)  J/D¥ 

CFREX32.*EE 

IF  (J.C. EC.  1)  GC  10  198 

If  (  (AES  (Cf  EjuL-CFREXJ/CFBEXL)  .GT.0.OQ5)  KB=3 
CFfct XL-CF  REX 
IF  (KB.  EC.  3J  GO  .TO  97 
198  KB=  3 

«  GO  1C  115 

150  JMAX=JKAXi1 

IF  (JKAX.LE.5Q)  GO  10  151 
.  JHAX=2b 

DO  116  J*1,  25 
JX=2*J  -1 

116  WL <1,J) =WL (1,JX) 

4L (1,2t)=4l  (  1,50) 

XECGE=¥EDGE+C? 

DX-C¥*2.0 
DESI=CESI*2.0 
CS1=DS1*2. 0 
DED»DEB*2.0 

i.  RI1 E (3,1 19)  08  AX,  Dt 

119  FORMAT  |/,3X,  '  NORMAL  GRID  POINTS  BECUCED  TO*  ,  14,  IX,/, 

*3X,'C (NORMAL)  = • , F 14. 6 ,/) 

GC  10  121 

151  rfRIlE  (2,1030) 

1030  FCRMA1  (5X,'1CINT  ADDED  AT  EDGE',/) 

¥££GEs¥£DGEfC? 

«L  (1  ,  JEAX)  =  4L  (1  ,JMAX-1) 

UEL*WI  (1,JMAX) 
xF  (UEI.IE.O.)  UEL= 1.0 
121  DC  113  J=  1 , JHAX 
113  i» (1.J) =21 (1,J)*UE/UEL 
GC  1C  £6 

115  FPE  =  -  (2. *4  (1  ,Jii  AX-  1) -0.5*4 (1.JMAX-2) -1. 5*4 (1, JHAX) ) /Dl 
CFREXE=CFBEX*FPE/FP 
IF  (NI.EC.  1)  GO  TO  117 

IF  <  (4  (1,JWAX-1)/9  ( 1 , JH  AX) ). IT. 0. 9995)  GOTO  150 
GO  TC  118 

117  IF (CFEEXE.Gt. 0.0005)  GOTO  150 

118  CONTINUE 

DO  120  J=1,JMAX 
■I (1,0) =4(1,0) 

120  CONTINUE 
H I  *  H 

IF  (CfREX.Gl.O.O)  GO  TO  130 

137  SS=-CFl/CCFD£  +  SI 

IF  (IK. EC. 0)  GC  TO  138 
•  XS=XL4CXDSl*  (SS-SL) 

PS*  (H4CPCS1*  (SS-SL)  )/DGB 
GO  TO  124 

138  PSl£=££/REEfi/OGR 
«R  1 1  fc  (3,1081)  PSIS 

1081  FORK A1  (10X, *  **  FLOW  SEPARATES  FOB  TuxS  STREAMLINE  AT  PS1  *',F7,3 
A,  IX, ' CEGREES' ,/) 

GO  TO  140 

124  »ii  11  E  (  3 , 1  C  80)  XS,PS 

1080  FiRiiAT  ( 10X, '  **  FLC2  SEPARATES  FOR  ISIS  STREAMLINE  AT  X  *',F7.4 
A,4X,'EKI  s  '  ,  F  7, 3, /) 

GC  TC  140 

130  if  (LP.Fc.O)  4RxT£  (3,1090)  (4  ( 1  , J) ,  J-  1 ,  J  MAX) 
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noon 


IF  (LK. FC.O)  GO  TO  131 
Y2C*  Y (2) /CGR 
Y3C=Y (3)/EGB 

DFES=£CEX*F  (1)-fDCPPH*F  (2) 

WHITE  (3,  1050)  Y  (  1)  ,  Y2C,Y3D,S,H,Ui,  OPUS  ,  CFiiLX 

IF(KIC.EQ.I)  WRITE  (3,  1052}  OP , DCPX , DCE PH , DCPX E, D2PPU, E , IT  EE, 

I  CF6EXF  ,  Y  (4}  ,Y(5) 

131  DCfCS- (CFSfi-CfL) /  (S-SL) 

CF I  =  CE E  EX 
S  L= S 
XL  =  Y  (1) 

PL  =  Y  (2) 

DXCSL=F  { 1 ) 

DPD£I  =  F  (2) 

IF (NHE.GE.HAXS)  GO  TO  140 
IF  (Y  (1) .GI.BL)  GO  TO  140 
IF(IK.EUO)  GO  TO  40 
GO  TO  550 

135  WRITE  (3 , 1 36) 

136  FORMAT  (/, 3X, 'LEST  PROFILE  FAILED  TO  CONVERGE  AFTER* , 

*1X,*40  ITEGATIQRS' ,/) 

GO  TC  137 

142  WR11E (3/143) 

143  FOE HAT (/ , 10X, 'STREAMLINE  TERMINATED  -  STEP  SIZE  REDUCED  TC  *, 

*  *  1CNEF  LIEU*,/) 

140  WRITE  (2/141)  NHP 

141  FORMAT (/, 3X, ' STREAMLINE  TERMINATING  AFTRH * , 14 ,  IX, *  ST AT10NS ' ,/) 
KBS*KES+1 

IF  (KfcS.LE.KEH)  GO  TO  0 

SICE 

END 


16.  Listing  of  Subroutine  BGEGM 

SUERCUTINF  BGEOM  (XX, S , DRCX, Gfl, DGX) 

C 

C  THIS  SUBROUTINE  COMPUTES  THE  GEOMETRIC  PROPERTIES 
C  FOR  A  S FU EBE-CG I VK-CY LINDER  CONFiGUR ATION 
C 

C  X  IS  AXIAL  POSITION 
C  B  IS  ECCY  RADIUS 

C  GEOMETRY  FOR  SPHERE-OGIVE-CYLINDER 
IF (XX. IE. C. 1442)  GO  TO  10 
If  (XX. GE. 5.  19615)  GO  TO  20 
R*-1 J.iSQRI (14. **2-  ( X X  —  5 .  19615) **2) 

DRCX- (5. 19615-XX) /  (R  +  1  3. ) 

D2RDx^a-(1.+CBDX**2)/ (R+lJ. ) 

GC  TC  30 

10  F  =  Sg51  (0. C 60011*4^  -  (XX-0,  1b6l43«) ***) 

DRCX* (C.  166 1 4 30-XX) /R 
D2RUX2*-  ( 1 ,  ■)  URD  X*  *2)  /R 
GC  TC  3C 
20  R*  1. 

Dl»EX«Q. 

D2RDX2*  0. C 
30  GM3ATAN  (DRCX) 
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n  n  r>  n 


DGX=C2BCX2/  (  1.4  DR  DX  *  *  2) 

RElUbN 

ENC 


SUBROUTINE  DUEOM  (X , S , DRDX , GM , DGX , TR AT) 

GECHElbX  ECR  AN  ELLIPSIOD  OF  REVOLUTION  WITH  THICKNESS 
R AllO  IbAl 

fi=Sgbl  (1RA1**2*  (1.0- (1.0- X) **2) ) 

DUUX=  (  1.0-X) *TR  AT**  2/R 
GM  =  ATAN  (DBOX) 

D2RDX2=-TRAl**2/&*  (1.  O-f(l.-X)  *CRDX/H) 

DGX  =  C2EDX2/  (1.0-)DRDX**2) 

RETURN 

ENC 


* 


17.  Listing  of  Subroutine  BLOTNR 


10 


11 


30 


SUBROUTINE  DLCT NR  |J?SI  ,  S  , X ,  F ,  DEXI ,  E X1H, BETA, DST) 
DIME  NS10N  X  (6) , F  (6)  , T  (4)  ,  P  (4) 

CCHMCN  /HA1L/D5TN,DX,UE 

COHHC  li  /OUTPI/B,  CP, DCPX,  DCPPH,  dcpxp, D2pph,gh, dgx 
COHflCN  /SCALt/HL,HLH,H 

CO  EMC b  /INTEC/XINT,PHIP,SAL#CAL,PI!,P5IO,BPEB 

IF (PS1.LE.PS10)  KL=0 
GO  10  (10,20.30) ,  KL 
T  1 1)  =0 . 0 
P(1)=0.0 

EQE=6EEB**2*RPER 

EX1L=0.0 

KL=  1 


T (4) =UE*H**2*MPEH 
P (4) =  P  £  1 
DISI  =  PSI-P  (1) 

DO  11  J=  1 , 2 
JK  =  J*1 

P(JK)  =E  (l)^DPSI /4. 0*2.0**  (J-1) 

CALL  SfHCAP  (P  (JK)  ,CPB,UED, DCP) 

IP  (JK. RE. 3)  GO  TO  11 

UCPS1=CCP 

UE3=UEE 

T(JK)  =EOB*Utb*SiN  (P  (JK)  )  **2 
X  (6)  =EXIL-)  CFS1*  (T(1)|4.*T(3)4T(4))/6.0 
EX1II=EX1L-)CPS1*  (T  ( 1 )  -f  4.  *T  ( 2 )  T  (3) )  / 1  2 . 0 
DEXI=X (fa) -EX1L 

BETA=-EXIH*DCPS1/  (T  (3)  *UE3**2) 

EX ILS X  (fa) 

P(1)  =  P  (4) 

T  (1)  =1  (4) 

UEL=UI 

IF  (PS1.GE.PH)  KL=2 
KE1UBN 

XH  =  Xl-)CXL*DSl/2.+DDXL*  (DST/2.)  **2/2.0 
X0=XL-|CXL*DSl/4.-|DDXL*  (DST/4.)  **2/2.0 
PH  =  mCEL*DSl/2.-fODPL*(UST/2.)  **2/2.0 
PC=PLFCPL*0Sl/4.+D0PL* (OSI/4.) **2/2.0 
CAU  IbVISDfXfi,  PQ,UEQ,DPX,DPP,  P,DP1P,DPP2) 
CALL  lb VI SC (XH,PU,UEH,DPX,OPP,P,DPXP,DPP2) 
HHM»MHL)/2.0  1 
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nonnnnnnnnonn 


lig=faI4  (h-HLj  /4.  0 

£UH=£XIL4Cai*(Ut’L*HL**2+4.*UE^*H(,**2  +  U£H*ilH**2)/U.O 

DXLS  =  EXL+CS1* CD  XL/2.0 

DPCS=CEL+£SI*CDfL/2.0 

BET A=-EXIH*  (CPX*DXCS4DPP*DPDS)  /  (  (UEi»*IIH)  **/*U£H) 

CtXI=X  (6) -EX1L 
EXIL  =  Y  (6) 

UEL=UE 
20  XL=X  ( 1 ) 

EL  =  X  (2) 

DXL  =  F  (1) 

DPl  =  f  (2) 

i)DXL  =  -SIN  (X  (J)  J  ^CGS  (Gil)  *1  (3) -COS  (X  (3))  *SXN  (Gil)  *DUX*F  (1) 
DDFI  =  CCS(X  (3))*F(3)/B-DBDX*F(1)  *SIN  (X{3))/R**2 
U  EL  =  U  £ 

KL*3 

FEIUuN 

EMC 


18.  Listing  of  Subroutine  FCN 

SUEbOUIINE  FCNJ  (N,S,X,PD) 

RETURN 

ENC 


JD2SCII1INE  FCM(W,S,X,F) 

DIMENSION  X  ( M)  ,  F (N) 

COfUCN  /hALL/CSIN,DX,UE 

cceacs  /outet/b ,ce,dcpx, dcpph,dcpxp,d2ppu,gji,dgx 

CC3ttCN  /5CAL£/hL,ttLU,H 

CALL  cGECF  (X  (1) ,fa,DREX,G3,DGX) 

CALL  IN VISC  (X  (1) , X ( 2) ,UE,DCPX, DCPEU, CP, DCPXP, D2PPU) 
CGf.=CCS  (G«) 

SG t  =  S I N  (GK) 

CIB  =  CCS  (X  (3) ) 

SIIi  =  SIN  {X  (3)  ) 

11H= i AN <X  (3) ) 

F  (1)  =ClK*CGfl 
F  (2) =STH/B 

F  (3)  =  0. 5*  jSIh*Cua*DCPX+CTH*DCPEH/F)  /UL**2-STH*SGM/F. 
F  (4)  =-X  (4)  ♦  (X  (4)+SG«)/ 

++0.  £♦  (STH*CGi1*DCPXP”CTll*D2EFH/R)  /UE**2 
•HO.  £  ♦CG«*CCEX*X  (4)  /  (CTH*UE**2j 

+4  0.E*ECEPn*  (SlH*CGa*DCPX-cTii*&CEPIi/«)  /(Uli**2*UE**2) 
F ( 5)  *X  (4) /  (n*CTH) 

IFjh.EC.5)  BETUBM 
li=S*Clk*EXE  (X  (S)  ) 
f  <c) =UE*H*** 
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BETUBN 

ENC 

C 

C 

C 

c 

c 

c 

c 

c 

c 

*  c 

c 

c 

c 


19.  listing  ot  Subroutine  COEFF 

SUEEOU1INE  CCiE'F  ( J  (1  AX  ,  Ail  , OH,  CH  ,  DM  ,  AS  ,  CS  ,  DS  ,  N  ,  ML,  KB) 

CCBMCN  /HALL/ES1N,DT,UE 
COMMON  /BLCl/CEIA.CEXI ,EXIH,BIPH 
ccanct.  /SCALE/ HL.HLH.H 
DIMENSION  h(2,JMAX)f  UL (2# JMAX) 

DIMENSION  Arf  (JMAX)  ,BM (JMAX) ,CM  (JMAX) ,DH  (JMAX) 

dimension  as (jm ax) ,os (jmaxj 
IF  (KB. EC* 2)  UO  10  00 
IF  (KB. EC*  3)  UC  10  80 
BE1A=C. 5 

C  SIBIL Ah  SOLUTION  COEFFICIENTS 

j  a  i= jm  ax-  i 

DC  20  J  =  2,  wMl 

AM (J  )  =  t .  5* (1 .0+0.5*DETA*h  (2, J) ) 

DM  (J)  =  1. 0-f  (BETA** 2)  *BETA*W  (1,J) 

Oil  (J)  =  1.0- Ail  (0) 

DM  (J)  =  C.5*  (DEI A** 2)  *  (BET  A*  (1. Q  +  N  (1,J)  **2)+W  (2,J)  *0.5* 

*(N  (1,Jf  1)  -«(  1,J-1))/DETA) 

AS  (J)=C.25*DtTA* (h (1.J-1) ) 

DS  (J) =C. 0 
20  CCMIUE 

C5=0. 5*CETA 
DS  (JMAX)  =  C.O 
EE1UHN 
C 

60  CONTINUE 
JM  1=JMAX-  1 
DO  0  5  J  —  2  r  J  il  1 

AS  (J)  =  C . S  +  C.  /5*D£1A*W  (2,  J) 

CM  (J)  =  1.0+ LtlA**2*li  (1,J)  *  (EIPIi+2.*EXIH/DEXI) 

C/i  (J)  =  1.0- AO  (J) 

DM  (JJ-C.5*  (  «  L  (1  ,J+  1)-2.*hl(1,J)+Wl  (1,0-1)) 

++0.  5*CElA**i*fclPH*  i»  (1,J)  ♦*2-«L  (1,  J)  **2+2.) 

+  +0.25*IET»*m  (2,  J)  *(K  (1#J+1)-H  (1,J-1)) 

++DE  lA**2*EXIii/DtXI*(*(  (1,  J)  **2+NL(1,J)**2) 

AS (J)  =C.  2 5* Cel  A* ( ta (1, J+1)-N ( 1, J-1)+WL  (1,J  +  1) -NL  (1.J-1) ) 
Co  ( J) =  CET  A*  (-0.25+LXIH/DEXI) *(hL(1,J)+UL(1,J-1)) 
to  5  CONTINUE 

iS-DETA*  (C .  ^5+  C.XI1I/CEXI) 

DS  (JMAX)  =  DEI  A*  (-0. 25+EXAIi/DEXI)  *  (2L  (1,  JHAX)  +  «L  (1.JHAX-1) ) 

f.  £  I U  h  N 
C 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

C  COEFFICIENTS  Fob  HALL*S  METHOD 
dO  CONTINUE 

U  (1,JMAX) =  UE 
J«1=JMAX-  1 
DC  90  J  =  2,Jrl  1 

AM  (J)=C.5/CY**2  +  0.25*b  (2,J)/CY 
BM  (J) =  1.0/CY**2  +  U (1 ,J) /DSTN 
CM  (J)3C.5/EY**2  -  W  (2,  J)  *0.  25/DY 

DM  (J)  =  0.5/C  Y**2  *{«!.(  1,  J+1)-  2.0  *irfL  <  1 ,  J) +WL  (  1,  J-1)  ) 

++(k  1,JJ  **2)  *0.  5/ DSTN 

++  (b  (1,JHAX) **2-WL(1,JHAX) **2)/  (2.0*DS1N) 

++0.25*fc  (2,J)/DY*(W  ) 

LS  (J) =CY*EL/  (2.  *DSIN*HLUJ  *  { »'  L  (1,J)-fhL  ( 1 ,  J  -  1 )  ) 

AS  ( J)  =  C.  25/C  Y*  (  N  1)-b  (1,J-1)+WL  (1,J  +  1)  ~UL  (1,  J-lj ) 

90  CCMINLE 

CS=CY»H/  (2«  *DSTN*HLU) 

CS(JMAX)=CY*UL/  (2.*DSTN*HLHJ  *  (hL  ( 1  , JM AX) +HL  ( 1  ,  JH4X- 1)  ) 
RETURN 
EMC 
C 
C 
C 
C 
C 
C 
C 

c 

c 

c 

c 

c 

c 

c 

c 


20.  listing  ot  Subroutine  INVERT 
susecutise  invert  (jm ax #a,l,c,d, as, cs,ds,«,kb) 

DIMENSION  A(JBAX),  b  (JflA  X)  ,C  (JHAX)  ,o  (JMAX) 
DIMENSION  AS  (JM AX)  ,  DS(JMAX) 

DIMENSION  U  (  2 , J  3  A  X ) 

Di-IEMSICN  E  (5C)  ,  EL('jO),  o  ( 30  ) 

COMMON  /HUL/CX,UY,Ue 
JMi=JM#X-2 
E  ( J M  A  X ) *0.0 
0  (JMAX)  =0.0 
EL  (J.4AX)  =  Ut 

IF  (KE.KE.  3)  EL  (JM  AX)  =1.0 
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DC  JK=  1 ,  JH2 
J  =  0!1AX-JK 

D£M  =  :»  (J)-l  (J)  *L  (u-)1)  +C3*  (C  (J)  *0  (J-)  1)  “  AS  ( J )  ) 

E  (J)  =  (A  (J)  -CS*  (C  (J)  *G  (J+1)-4;>  (J))  )  /DEM 
(j  (J)  =  (C  (J)  *0  <JF  1)  -AS  (J) )  /DEM 

LL  (J)=  (C  (0)  +  (C(J)  *C  (J+1)  -AS  (J)  )  ♦DS  (o)+C  (J)  *EL  (J+ 1 )  )  /DEM 
20  CONTINUE 

DO  10  J  =  2,J11A1 

f  U  (1,0)  =E  (J)  *k  (1,J-1)+G(J)  *W  (2.J-1)  +EL(J) 

M(2,J)=U(2,J-1)-CS*(M(1,J)+2(1,J-1)j+D3(J) 

10  CCMHMCE 

t  ME11JSN 

lUL 

c 

c 

c 

c 

c 

c 

c 

c 

c 


21.  Hating  of  Function  KliUMGE 
FUNCUCN  KFUNOE  (X  ,  P  #  L  ,  H,  »  ,  H.1AX,  MB) 

DIKE  AS  IC  N  P 11 X  <  1  2)  , SAVE! (12)  ,  X  1  (12)  ,  X2  ( 1 2>  ,  XKP  1 12)  ,  PKP  ( 12)  ,X(M) 
*F(N) 

DAI  A  E!,LOCF,fcPa/0,0,5. £-4/ 

M=K+  1 

CO  10  (5,45,65.65) ,rt 
5  IF  (ICCf.GI.O)  GO  TO  25 
IF  (ilfc.Eu.  1)  00  IC  205 
IF  (AES (H) .OF. ADS  (UMAX))  U  =  HBAX 
DO  15  J=1,N 
XKt  ( J) -  X  ( J ) 

15  FKE (J) * F  (J) 

XO  =  X 

25  DC  35  J  =  1  ,  N 
SAVE  i  (J) =  X  (0) 

Pill  (J)  -F  (J) 

35  X  (J)  =  S  AVEX  (0)30.5*H*F(J) 

X*x-t0. 5*f( 

KBUNGE=  1 
hETUNN 

45  DC  55  J  =  1  ,  N 

Pill  (J)  *FIII  (0>+2. 0*F  (J) 

55  X (J) =SAVEX (J) +0.5*H*F (J) 

uElUftN 

65  DC  75  J=1,N 

Pill  (J)  =  fiii  (J)  *2.0*F  (J) 

75  X (0) =SA VEX  (J) jb  +  F  (J) 

X  =  X-)0 . 5  *H 
hElUMN 

«5  DC  95  J  =  1  ,  M 

95  X  (J)  =SAVEX  (C)-l(Phl  (J)+F(J,  )  *H/6  •  0 
IF  (MF.Efi.  1)  00  1C  165 
IF  (LCCt-  1)  105,  125,  145 

105  DC  115  J=1,K 

X 2 (J)  =X  (J) 


125 


n  n  n  n  n 


e  (0) =Fht  (J) 

115  Y (0) =  YKE (J) 

X  =  J<C 
ii  =  h/ 2 . 

a=  1 

ICCF  = 1 
GO  1C  25 

125  DC  135  0=1, il 
135  Y  1  (0) =  Y (J) 

LCCP=2 
M  =  0 

KElUfcK 

145  DC  155  0=1, N 

*F  ( ALS (Y  (J) j .LX. 5. £-04)  GO  TO  155 
EB=  (Y  (0)-Y2(J))/Y  (J) 
it  ( A  £  S  (E  R ) -  EPS)  155,  155,  175 
155  CCMIMJE 
165  H=2. *H 
MB  =  0 
LCCE  =  0 
KHUNGE*0 

a=o 

8E.1UKN 

175  DC  185  0=1,41 
Y (J)  =  YKF  (J) 

F (C) =  FKE  (0) 

185  Y2(0)=11(0) 

X=XO 

H=H/2. 

if  (AtS  (»)  .11. 1.0-10)  GO  XO  195 

LCC  F®  1 

a=i 

GO  10  25 
195  K6UNGE=2 
a=o 

LOCF  =  0 
RETURN 

205  DO  215  0=1, il 
i (J) =  YKP  (0) 

215  i  (Gj  =  F  K  £  (J) 

X=i  C 
n-ii/ i • 

GO  10  25 
ENt 


22.  listing  o £  Subroutiue  milneg 

GUt ROUTINE  MiLNEG  (Y,F,Pca,.NDE,DG,S) 
DIMENSION  Y (6) ,r (o)  ,PC;i(o,4,  2) 

5  DC  1C  0=1,  NCE 

Y  (J)  =FCK(0,4, 1)  +-».*DE/3.  *  i2.*PCK(J,  1,2) 
++  * . *  C  C  K  (0,3,2)) 

10  C  C  n  1 1 M  E 

CALI  E L  N  (  N Ck  ,  8, Y , F) 

DC  30  0*1,  NCk 


“PCM (0,2,2) 
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n  r,  n  n  r.  n  n  n  <\  n  r>  c.  n  n  no  n 


X  (J>  =  t  C  K  (J.2, 1)  (-OS/J.  ♦  (t  (J)  iU.  *PC.1  (J,  1,2)  -f  PCM  (J,  2,2)) 
30  v.C  NT  I  N  L  E 

CALL  tCi  (NEt,S, V,F) 

DC  bO  J=  1,  NCE 
DC  6C  JJ=1,J 
JK=-S-J  J 
JK2=JK-  1 
DC  7  C  JL-  1 , 2 

70  PCfl (J,  JK,  Jl) =  tCfi  (J,JK2,JL) 

60  CONTINUE 

PcM  (J,  1 ,  1)  =X  (J) 

50  pcm  |J,  1,2) -t  i j) 

S=S(DS 

RETURN 

LNC 


23.  Listing  o£  Subroutine  PRESS 
SUE ROUTINE  PRESS 

CUBUCN  /PCM  A/P  Hi  (40)  ,X  (40)  ,  CP  (40,40)  ,CPX  (2,40) 
CCCBCN  /CtST  A/N AS,NCS 


NAS  =  N  U  H  L  h  R  C t  AXIAL  PRESSURE  STATIONS 

NCS  =  NUtUER  CE  CiRCUBLERENTxAL  PRESSURE  STATIONS 

Du  R  =  AC  C  S (-  1. ) /I  oO. 

REAL  (1,5)  NC 3 , NAS 
5  FOPflAl  (12,  IX, 12) 

DC  10  J  = 1 , Nc  S 

REAL  (1,  15)  tlil(J) 

PHI  (J)  =PI11  (J)  *UGR 
15  ECIiflAT  (FlO.b) 

10  CONTINUE 

DO  xO  J  =  1 , N A S 

HEAD  (1,15)  X  (J) 

20  CONTINUE 

DC  30  J  =  1  ,  N  A  S 

RE  AD ( 1 , 35)  (CP(J,K),  K= 1 , NCS) 

35  FCRBAT  (5  (  1X,HO.'j)  ) 

30  CONTINUE 


0PX(1,J)  AND  C P X  ( 2 , J )  ARE  Tliii  AXIAL  DEolVATlVES  CF  THE  PRESSURE 
CuEr  E  ICIEST  AT  IHl  INTERFACE  AND  RCDX  END.  EACH  CIRCUMFERENTIAL 
PLANE  ML ST  HAVE  A  SET  AS  INPUT  TO  TuE  WU  ADR AIIC  SLPINE  ROUTINE. 

DC  UO  J  -  1 , NC S 

TEAL  (1,  16)  ctX  (  1,  J)  ,CPX  (2,J) 

1 6  FCRK  AT  (F 10. b,  IX, t  10.5) 

40  CONTINUE 
RETURN 
ENE 
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n  n 


24.  Listing  o t  Subroutine  SPHCAP 

SUDBCUIINE  SPHCAP  (P,CP,U£,DCFSI) 

DIME  NSICN  COE  (10) 

DAI  A  CCE/0. 666823,  1.  1 30£>72#  -0.1343b0,  0.  131321,  -0.009964, 
*  0.014246,  C. 016485,  -0.003913,  0.014531,  -0.627853/ 

CCCflCA  /IA7fC/X£NT,PfiZP,5AL,CAL,Pn,P5IO,RPE£ 

If  (f • E£. 0 , C)  M=0 
GC  1C  (210,200)  ,» 

CALL  IMVISD (X1NT,PH1P,U£, DCPX2 , DcP EH2,CP2, DCPXP, D2PPH) 

CALL  EGECfi  (1INT,B,CR,G.'1,  COX) 

DF3IEX-  (SAL*S1N  (GM)  ♦COS (PHIP)  -CALICOS  (Gtl)  )  ♦DOX/SIN  (PH) 
CCEDSI=G. C 
C  E=  0 . 0 
DO  10  J=  1 , 9 

cp=cfgcc£ ( Jj ♦cos(J^Prt) 

10  CCEDS  J=CCEESl-J*CO£  (J)  *SIN  (J*PM) 

CP=CP-fCCE  (10) 

DCEX=CCEDSI*CPSIDX 
2 8 1 1 £  ( J,2C)  CP,CP2,DCPX,DCPX2 

20  FCB.VA1  (/,4X, 'CONTINUITY  Of  PHES3UBE  AND  DERIVATIVE  ACROSS', 
A*  IN1EEE ACE:  • , / , 2  (5X , F 1 2 . 5) , / , 2  1 5X , F 1 2 . 5)  , /) 

.1=2 

IF  (ALS  (CP-CE2)  .GT. 0.00 3. Oh.  ADS  (DCPX-DCPX2)  .GX.  4.0)  11=  1 
IF(f!.£C.1)  NhlTE  (3,30) 

30  FOBflAI  (4X,  '**•**  QUADRATIC  SPLINE  EXTENDED  TO  INCLUDE', 

A*  NCSE  BEG IC  N  *****',/) 

IF(fi.EC,2)  B  ETUhN 

PH=EK/2.0 

CFE=LCFX2/EFSIDX 

Y=CE2/2.04C.5-CP?*FH/8.0 

YP  1=6.  C*  ( Y-  1 . 0) /PN**2 

YP2  =  b.C*  (Y-CE2+CPP*Pt1/2.  0)  /PS* *2 

fiEIUEN 

210  IF(E.Gl.Ffi)  00  Io  220 
CE=1.C-)YP1*t*  *2/2.0 
DC  ES 1=  Y  E  1  *  P 
OG  10  230 

220  CF  =  CE2-ECPE*  (f-Pf!)  -f  YP2*  (P-P«)  **2/2.0 
DCESI  =  CEP-)YE2*  (P-PS) 

230  UE=5CB1  (1. 0-CE) 

RETURN 
200  CF=0.C 

DCFS  1=  C. 0 
DO  250  J=  1 , 9 

ECESi=CCPSI-J*CUE (J) *SIN (J*P) 

250  CP  =  CE-)CCE  (JJ  *003  (J*P) 

Ct  =  C  E-f  CCE  (10) 

UE  =  SCB1  (1.0-CE) 

REIURN 

END 

C 

C 

c 

c 

c 

c 
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25.  listing  of  subroutine  (llEPTS 
SUESCllINE  (1IEP1S (AA, BB,LC,D, YS,N) 

OIPikSltN  C  (40)  ,0  (40)  ,  Kk (40) ,DB  (40) ,CC  (40) ,D  (40) ,  YM  (4  0) 

hH 1  =  N-  1 

NMk=N-2 

t  (1)  =  -CC  (1)/LE(1) 

U(1)=C(1)/t1:(1) 

DC  5  K=2,MB1 

ALEhA-AA  (K)  *(,  (K-l)-)liB(K) 

W  (K)  =-(.(.  (K)/ ALPHA 

5  U  (K)  *  (C  (K)-AA  (K)  *U  (K-1)j  /  ALPHA 
Yfl  (NB1)  =U  |NH1) 

DO  10  K=*1,AB2 
JK«A!l2-(  1-K 

10  YB  (JK)=C(JK)  *YB  (JK+1)+0(JK) 

EtEIUBN 

EMC 


26.  listiuij  ot  Subroutine  IMViSO 

SUEBCUl JNE  Ii«VlSC(XX,PPH,Ufc,DLPX,Da  PH,CPC,DCPXP,  D2PPH) 
CCBBCA  /ECU  A/PUI  (40)  ,X{40)  ,CP  (40,40)  ,CPX  (2,40) 

ccnncN  /cisi a/n as, nos 

DIRELSICN  CL LX (40) , ££LP (40) , DP  (40) ,CP1  (40,40) ,CP2  (40,4  0) 
DldeNSlCN  A  (40)  ,0(40)  ,  C  (  40)  ,  U  (40)  ,  Yil  (40)  ,  CPT  (2,40) 

LlREKS  ICN  ALP  (40)  ,  BET  A  (4  0)  , GAB  (40) , DEL  (40) 

DA  I A  t./ 0/ 

IF  (M. NE.O)  OC  TO  100 


DO  10  J=2,N*S 
10  DELX  (J)-X  <J)  -X(J-  1) 
NCR  1  =  NCS- 1 
ULH2=MtS-2 
MF1*NAS-1 
NAB2-MAS-2 
A  (1)  =C.C 
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B  (1)  -2./0EL1  UH  (2.4DELX  (  J)  /DELX  (2) )  /  (DELX  (3)  +CELX  (2)  ) 

C  (1)  =  C  E  LX  (2)  /  DEi.X  (3)  /  (DtLX(2)+DELX  (J)  ) 

DO  20  £=2, NAH2 

A  ( N }  - C E  1  X  lSi1)/DELX  (N)  /(DtLX  (Nil)  +  DELX  ( N)  ) 

B  (N)  =  (2.4  CEtX  (N )  /DELX  (N+l)  )/  (DtLX  (N)  40ELX  <N+  1 )  ) 

44  ( 2,. 4 C E L X  |N42)/DELX(N4l))/(DELX  (N42)  4DELX  (N4  1)  ) 

C  (N)  =  CELX  (Nt  1 )  /DELI  ( N  4  2)  /  (DELX  (N4  1)  fDtLX  (Ni2)  ) 

20  CCS1IN0E 

A  (NAK1) -DELX  (NAS) /DELX (N AM  1 ) / (DELX  (NAS) 4  DELX (NAH1) ) 

Q  (NAM  1)  =  (2. 4 DELX  ( « Aft  1 )  /DtLX  (NAS)  )  /  (L)£LX  (NAS)  4  DELX  (NAM  1)  ) 

442. /DELX (NAS) 

C(SA«1)=0.0 

C 

C 

DC  30  N- 1 , NLS 

D(1)=C£  1 1  #  M)  *2./DELX(2)+CE  (2,  N)  *(  (2.  4  DELX  (3)  /DELX  (2)  )  /  (DELX  (3) 
44-DEL  X  (2) )  4CELX  (  2)  /  DELX  (3)  /(DELX(2)4DEEX(3)  ))  4  CPX(1,N)/2.0 
DO  32  J-2.NAM2 

D  (J)  =  C£  (J,N)  *  (A  (0)4(2. 4DLLX  (J) /DELX  (J4l)  )/  (DELX  (J)  4DELX  (J41)  ) ) 
44CE  (J4  1  ,N)  ♦  (C(d)  4(2.4DELX(J42)  /DELX  (J4 1 ) )  /  (DELX  (J42)  4DELX  (04»)  >) 
32  CCN1INUE 

D (kAHl)  =CE (N AN1,N) * (DELX (NAS) /DELX  (NAM1)/  (DELX  (NAS) 4DELX  (NAM  1) ) 
44  (2.4DELX  (NANI)  /DELX  (NAS)  )/  (CELX  (  NAM  1)  fCELX  (NAS)  )  ) 

44C £  (NAS, N) *2. /DELX  (NAS)  -  CPX(2,N)/2.0 
CALL  H I  DPI S  (A,B,C,C,YH,N  AS) 

DO  34  J-2f  NAM  1 
11-CLLX  (J41) 4DELX  (J) 

12  =  (XM  (J)-CE  (J,N)  )/DELX(J4l) 

13"  (YR  (J-1)-CP(J,N)  )  /  DELX  (J) 

CP  1  (J ,  N)  *2. /II*  (DELX(J)  *1 2-* DELX  (J4I)  *TJ) 

CP2(J,A)=>«./1 1*  (12413) 

34  CCNI1NUE 

CPI  (1,N)=CEX  <1,N) 

CP2(1,N)*e.*  (YM  (1)-CP(1,N)-CPX  (1,N)*JELX  (2)/2.)/DELX(2)  **2 
CPI (NAS, N) =CPX (2«  N) 

CP2  (AAS,N)*S.*(YH  (NAM  1) -CP  (NAS.N) 4CPX  (2,S) ♦DELX (NAS)/2. )/ 

ADELX  (NAS) **2 
30  CCkllkUE 
C 
C 

DC  35  J=1,NAM1 

35  DELX  (J)  =  (X  (Oil)  4X(J)) /2. 

C 

DC  40  J-2.NCS 

40  DELE  (J)=PBL  (J)-PHL  (J-1) 

C 

Alf  (1)  *0.0 

BE1A  (1) *2./UtLP  (2) 

(JAM  ( 1)»  (2.4CELP  (3)/DELP(2))/(DELP(3)  4DELP(2)  ) 

DEL  (1) =  LELP  (2) /DLLP (3) / (DLLP (2) 4OELP (3)  ) 

A(1)*Alf  (1) 

D  ( 1)  =EEIA  <  I)  4OAM  ( 1) 

C  (1) =EEL  ( 1) 

DO  45  J*2,NlM2 

ALP  (J) *CEIP  (J41)/DELP  (J)  /  (OELP  (J4I)  4DtLP  (J)  ) 

BETA  (J)  =  (2 .4  CELP  (J)/DELP  (J4l) )/  (DELE  (J)  4DELP  (Jfl) } 

SA3  (J)  =  (2.4LELP(J42)/DELP  (J4I)  )/  (DELP  (J42)  fDiLP(J4  1)  ) 

DFi  (J)-CELP  (J4l)/DELP  (J42)  /  (DuLP  (J 4  1)  40ELP  (J 4 2) ) 

A  (J)  =  AL  f  (J) 
u  (J)  =  E E 1 A  (J)  40 AH  (J) 

C (J) *CEL  (J) 
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45  CCAUALfc 

ALP  (NCR  1)  -ZLlL  (NC3J/DELI  ( ACfll )  /  (DELP  (DCS)  +  DELP  (NCfl  1 ) ) 

BETA  (ACfll)  =  (2.+DLLP  (NCM)/DELP (NCS) ) / (DELP (NCfll)  +DELP (NCS) ) 

GAB  (ACM)  =2. /DELP  (KCS) 

DEL  (ACK  1)  *C.O 
A  (ACM  1 )  =  A L F  (NCH  1) 

MAC  Ml)  =DEIA  (NCdl)  +GAfl  (NCdl) 

C (ACfll) =  DEl  ( ACfl  1) 
c 
c 

DO  50  JM, ACfll 

50  DP  (J)  =  (flU  (J  +  1)  +  PH1  (J)  )/2.0 
fl=1 
C 

c 

c 

c 

c 

100  DO  105  J~  1 , A Afl 1 

IF  (XX.LE.EELX(J))  GO  TO  110 
105  CCA1IACE 

110  IF (XX. Cl. CEI*  (NAfll) )  J=NAS 
DEX=XX-X  (J) 

DC  115  K=  1 , ACS 

CFI  (1,K)=CP(J,K)+CP1  (J ,K)  *DiiXfCP2  (J.K)  *D EX**2/2. 

115  CFI  (2  #K)-CF  1  (J,K)+CF2  (J,K)*DEX 
C 

C  SPLINE  FIT  DC1H  CP  AND  E(CP)/D(X) 

C 

DO  1J0  LM, ACfll 
IF  (FEH.LE.EF  (L)  )  GO  TO  140 
1JO  CONTINUE 

140  IF  (PEL. GT. CP  (ACfll) J  L*NCS 
DFKI  =  I EH-EH1  (L) 

C 

DO  12C  A«  1,2 

D  (1)  =  CE1  (N  ,  1 )  *UET A  ( 1)  +CPT  (N.2) *  (GAS  (1)  f  DEL  ( 1)  ) 

DC  125  K*2,  ACB2 

125  0  (K)  =CET  (A,  A)  *  (ALP  (K)  ■(■  BE  T  A  (K)  )  +CFT  (A,  K+1)  *  (GAfl  (X)  -f  DEL  (K)  ) 
MNCH1)  =CET  (N, ACfl  1)  *  (ALP  (NCN1)  (-BETA  (NCfl  1)  )  +CPT  (A,  NCS)  *G  Afl  (NCM1 ) 
CALL  flICPIS  (A,d,C,E, Yfl.NCS) 

If  (I. Ft.  1)  GC  TO  12b 
IF  (L.EC.NCS)  GO  TO  127 
1 1=CELP  (1  + 1) +CELP  (L) 

T2=  ( Yfl  (L)  -CPI  (N  ,L)  )  /DELP  (L+1) 

13  =  ( Y  K  (l-l)-CPT  (N,L) J/DELP  (L) 

TJ=2./T1* (CtLE(L) *T2- DELP (L+1) *T3) 

1K=8./11*  (12+13) 

GC  TC  128 

126  T  J  =  0 .  C 

Tri  =  8.  * (Yfl (1)  -CPT ( N  ,  L) ) /vblP  (2)  **2 
GC  TC  Ud 

127  TJ=C . C 

TK=8.»  (Yfl  (ACfll)  -CPT  (N.L)  )/DELP  (NCS)  **2 

128  CPI (N, 1) *Ctl  (tl.L) +lJ*DPUl+TK*CPHl**2/2. 

GC  IC  ( 122.  1i 4)  .  N 

122  IF (Prk.Cl.O. )  GO  TO  123 
DCP  c  H*  G  .  0 

L>2FEi)  =  1K 
GC  1C  120 

123  DC  FI H=TJ+TK*CIHI 


n  n  n  n  on  non  non 


02PPh=TK 
GC  IC  120 

124  DCEXP=TJ+TK*CPHi 

120  CCN1INUE 

CPC=CFT  (1,1) 
LCEX=CE1  (2,1) 
UE=SCE1 (1.-CPC) 
RETURN 
ENC 


SUchOOTINE  INVISD(X,PHI,CP,UE,UEX,OEP) 

PCIENTIAI  SOIOTICN  FCfi  ELLIPSOID  Of  BEVOLUTION 

CCBMCN  /SlAu/XO, XC,3AL,CAL,TRAT,SC,OBE.SBE,CGO,AO,BO, VGA  AD 
DATA  B/0/ 

IF  (fl.EC.  1)  GO  TO  10 
11=2. 0*CAI/ (2.0-AC) 

12*2. C*SAI/ (2.0-BO) 

10  CALL  EGECH  (X,S, EHDX,GH,DGX,TBAT) 

APSC  =  (A- 1.0)  **2-)  (fi/TRAT**2)  **2 

15=  (1 1*B/TSA1**2+T2*  (X-1. 0)  *CCS  (PEI)  ) 

UE=SCRT (T5**2/APSU+(T2*SIN  (PHI))**2) 

CP= 1 . 0- UE**2 

U  EX*  (15*  (1 1*CRDX/TRAI**2+I2*CCS (PHI) ) /APSu>T5**2* (  (X-  1 . 0)  f 
A  H*DRDX/TBA1**2/T8AT**2) /APSC**2)/0b 
UEF=  (-15*  (T2*  (X-1. 0) *SIN  (PHI) ) /APS()+T2**2*SIN  (PHI) *COS (PHI)  )  /UE 
If(B.EC.I)  RETURN 

DUSCX*  (1 1  *  EBCX/TR AT**2+T2) /SORT (APSw) 

VGRAO=EUSCX*CGO/T2*XO 

a-i 

REIURN 

END 


27.  Listing  ot  Subroutine  STAGE 

SUBROUTINE  STAGE ( ALP, XO, XNOSK, HPEP) 

THIS  SUBROUTINE  COMPUTES  THE  NEWTONIAN  STAGNATION  POINT  FOR 
A  GIVEN  EffECTIVE  ANGLE  OF  ATTACK 
XC=XNCSE+RFEE* ( 1.  0-COS (ALP) ) 

RETURN 

END 

SUE RCU TINE  STAGN (ALP,XO, TO, TH AT, AC , BO) 

THIS  60U1 1 NE  LOCATES  THE  STAGNATION  POINT  CN 
AN  ELLIPSICD  CF  REVOLUTION 

E*SCR1  (1. 0-TRAT**2) 

ES£=E*£ 

£03=  ESC 

AO* 2.0* (1 .0-ESQ) * (0.5*ALOG((1.+E)/  (1. -E))-E) /EyB 
DC*  1 .  C/ E SC**  (1.0-ESo)  *  ALOG  ( ( 1 .  +  E)  /  1 1.  - E)  ) /  (2.  0*E0B) 

E*  (2.C-AO) *1  AN (ALP) /  (2.0-u0) 

XO*P*T RAT **2/SU hi  (1.0+  (P*1RAT) **2) 

XO=1.0-SGB1  (1.0- (XO/TRAT) **2) 

RETURN 

ENC 
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28.  Listing  01  Input  Data  tor  Sphere- 
c^i ve-Cylindcr  Configuration 

0.060811  0.  1053328  17.50  0.1442 

0.C01U  2.  0  65.0 

000  000  0C1  003  300  000  COO  000  002  COO  000 
10  34 

,  0.0 

11.25 
33.  75 
52.  25 

*  76.75 

101. 25 

123.75 

146.25 

168.75 
180.00 
0. 1442 
0.  15690 
0.  16 

.  0.17 

0.  18 
0.20 
0.21630 
0.26620 
0.32870 
0. 40964 
0.50 
0.  60 
0.70 
0.80 
0.90 
1.00 
1.  100 
1.20 
1.50 

1.75 

2.0 

2.5 
3.0 

3.5 
4.  0 

4.5 
4.61 

l  5. 66 

6.97 
8.  29 

*  9.61 
12.  24 
14.87 
17.  5 


0 • 98200  COO 

G.058COOOO 

0.  16400000 

-0.4890O00 

-1.3840000 

-  1. 6bH 

-1.56c 

-1.050 

-1. oJU 

-1. 598 

0. 9716 

C .  6  6  9  o 

0.2320 

-0. 7503 

-1.5789 

-  1. 694  7 

-  1. 6499 

-1.  1143 

-0. 6985 

-0. 6465 

C .  9  7  0  3 

0.6933 

C.24C8 

-0. 7800 

-1.610 

-  1. 90  10 

-  1 . c3o0 

-1.09 

-0.695 

-0. 639 

0. 9605 

C .  8  5  8  9 

0.2570 

-0.6450 

-1.6620 

-1.31o 

-  1. 604 

-1.012 

-0.6/2 

-0.  6  1  4 

0.  905 

C  .  900 

0.26  1 

-0.6540 

-1.664 
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-1.917 

-  1.  580 

-0.9bS 

-0.653 

0.9701 

C.9CJ0 

0.2403 

-0.8330 

-1.891 

-1. 5440 

-0.904 

-0.619 

0.9770 

C . 9047 

0.2193 

-0.8100 

-1.8723 

-1.5214 

-0.6717 

-0.5950 

0.9890 

C  .9040 

0.2242 

-0.7777 

-1. 6328 

-1.4710 

-0. 8134 

-0. 5300 

C. 9862 

C . 9027 

0.2347 

-0.7498 

-1. 76b2 

-1. 4304 

-0. 784 

-0.4700 

0.9811 

C.8995 

0.2460 

-0.7169 

-1. 7303 

-1. 3819 

-0. 7492 

-0.4200 

0.975 

C  .894 

0.256 

-0.68 

-1.7 

-1.  22» 

-0. 725 

-0. 382 

0.967 

C.886 

0.262 

-0.656 

- 1. 666 

-1.262 

-0. 694 

-0.354 

0.  958 

C.882 

0.267 

-0 . b39 

-1. 636 

-1.235 

-0.667 

-0. 335 

0.95 

C.874 

0.269 

-0. 626 

-1.608 

-1.  190 

-0.636 

-0. 321 

0.942 

C.867 

0.  271 

-0.615 

-1.583 

-1.  147 

-0.608 

-0. 314 

0.935 

C.859 

0.272 

-0.606 

-1.55S 

-1.  105 

-0.580 

-0. 309 

0.  92b 

C.85  1 

0.271 

-0.597 

-1.  536 

-1.  C68 

-0.552 

-0. 305 

0.918 

C.842 

0.270 

-0.592 

-1.516 

-1.C33 

-0.525 

-0. 300 

0.894 

C.81S 

0.259 

-0.580 

-1. 467 

-C. S50 

-0.450 

-0. z8  6 

0.872 

0.794 

0.243 

-0.575 

-1.435 

-0.693 

-0.391 

-0. 274 

0.849 

C.770 

0.223 

-0.572 

-1.  412 

-0. 634 

-0. 337 

-0.262 

0.798 

C.72Q 

0.  176 

-0.574 

-1. 384 

-0. 7*1 

-0, 253 

-0. 245 

0.741 

C.666 

0.129 

-0.582 

-1. 347 

-0.624 

-0. 194 

-0. 205 

0.  604 

C.613 

0.080 

-0.595 

-1. 284 

-0. 54a 

-0. 154 

-0.  178 

0.627 

C.563 

0.033 

-0.610 

-1.  199 

-0.  477 

-0.  162 

-0.  156 

0.583 

C.518 

-0.0  11 

-0. 627 

-1.  105 

-Q. 415 

-0.215 

-0.  147 

0.5570 

C • 510b 

-0.0196 

-O.bJOS 

-1.083 

-0. 3964 

-0.  23 

-0. 150 

0.5441 

C . 4730 

-0.0528 

-0.6705 

-0.860 

-0.  39 

-0.425 

-0. 4477 

0. 5380 

0.4639 

-0.0790 

-0.7184 

-0. 7187 

-0. 4 622 

-0.7553 

-0. 8456 

0.531 

0.4526 

-0. lOOo 

-0.7546 

-0.609 

-0.5472 

-0. 858b 

-1. 0870 

0.525 

C.4465 

-0. 1148 

-0.7658 

-0.  84 

-0.635 

-0. 858 

-1. 189b 

0. 5144 

C  .4228 

-0. 1205 

-0.7538 

-1.067 

-0. 64o  3 

-0.6351 

-0. 7362 

0. 5080 

0.4226 

-0.074b 

-0.6380 

-0.54 

-0.  22 

-0.  38 

-0. 4S3b 

0.5012 

C.43S4 

-0.0131 

-0. 5321 

-0. 415 

-C. 280 

-0. 340 

-0. 3949 

■1.480 

-C. 01225 

i.  60 

C.CCC74 

-0.595 
-1.6310 
-0. 5  o  0 
-1.6163 
-0. 5380 
-1.5840 
-0. 4650 
-1. 5430 
-0.  4300 
-1.  4910 
-0. 3620 

-1.432  * 

-0. 345 
-1. 398 

-0.317  • 

-1.374 

-0. 300 

-1. 355 

-0. 284 

-1.34 

-0. 277 

-1.326 

-0.274 

-1.315 

-0.272 

-1.306 

-0. 2b9 

-1.287 

-0. 2b4 

-1.278 

-0. 256 

-1.272 

-0. 254 

-1.264 

-0.245 

-1.264 

-0. 224 

-1.277 

-0. 177 

-1.290 

-0. 148 

-1.299 

-0. 137 

-1.3070 

-0. 140 

-1. 3163 

-0.  450 

-1.2856  j 

-0.85 

-1.  1897 

-1.117 

-1.2440  * 

-1. 240 

-1.4508 

-0. 750 

-1,0941 

-0.509 

-0. 8233 

-0. 402 
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3.40 

C.C2270 

-16.67 

• C . 04  60 

-15.0 

C.C804 

-2.  20 

C.C1683 

15.0 

-C. 01633 

32.  25 

-C. 0G4d J 

47.5 

C.C130 

53.  25 

C.C1633 

50.  0 
75.0 


l 


m 
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